% Create a format configuration object
cfgDMG = wlanDMGConfig;
mcs = "9"; % MCS specified as a string scalar or string vector
cfgDMG.PSDULength = 4096; % PSDULength in bytes

% Get sampling rate and specify carrier frequency
fs = wlanSampleRate(cfgDMG);
fc = 60e9;

% Create a TGay channel object
tgayChan = wlanTGayChannel;
tgayChan.SampleRate                = fs;
tgayChan.CarrierFrequency          = fc;
tgayChan.Environment               = 'Open area hotspot';
tgayChan.TransmitArray.Size        = [4 4];
tgayChan.TransmitArrayPosition     = [0; 0; 6];     % Meters
tgayChan.TransmitArrayOrientation  = [0; 270; 0];   % Degrees
tgayChan.ReceiveArray.Size         = [4 4];
tgayChan.ReceiveArrayPosition      = [6; 6; 1.5];   % Meters
tgayChan.ReceiveArrayOrientation   = [90; 0; 0];    % Degrees
tgayChan.BeamformingMethod         = 'Maximum power ray';
tgayChan.NormalizeImpulseResponses = true;
tgayChan.NormalizeChannelOutputs   = false;

showEnvironment(tgayChan);

ppm = 20;                 % Clock accuracy to drive the CFO (ppm)
freqOffset = ppm*1e-6*fc; % Carrier frequency offset (Hz)
delay = 500;              % Sample to delay the waveform
zeroPadding = 100;        % Add trailing zeros to allow for channel delay

snrRanges = {-2.0:0.5:0.5, ...  % MCS 1
             0.0:1.0:5.0, ...   % MCS 2
             1.0:1.2:7.0, ...   % MCS 3
             2.0:1.2:8.0, ...   % MCS 4
             2.0:1.3:8.5, ...   % MCS 5
             2.5:1.3:9.0, ...   % MCS 6
             4.0:1.3:10.5, ...  % MCS 7
             5.0:1.5:12.5,...   % MCS 8
             5.5:1.5:13.0, ...  % MCS 9
             7.0:1.5:14.5, ...  % MCS 9.1
             8.0:1.8:17.0, ...  % MCS 10
             10.0:2.0:20.0, ... % MCS 11
             12.0:2.0:22.0, ... % MCS 12
             12.0:2.0:22.0, ... % MCS 12.1
             14.0:2.0:24.0, ... % MCS 12.2
             16.0:2.5:28.5, ... % MCS 12.3
             17.0:2.5:29.5, ... % MCS 12.4
             17.0:2.5:29.5, ... % MCS 12.5
             20.0:2.5:32.5};    % MCS 12.6

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

numSNR = numel(snrRanges{1}); % Number of SNR points
if ~isstring(mcs)
    error('MCS must be specified as a string scalar or string vector');
end
numMCS = numel(mcs);          % Number of MCS
packetErrorRate = zeros(numMCS,numSNR);
Ngi = 64; % Fixed GI length defined in the standard (20.6.3.2.5)
validMCS = string(sort([1:12 9.1 12.1:0.1:12.6]));

for imcs = 1:numMCS
    cfgDMG.MCS = mcs(imcs);
    if ~strcmp(phyType(cfgDMG),'SC')
        error('This example only supports DMG SC PHY simulation');
    end
    ind = wlanFieldIndices(cfgDMG);
    snr = snrRanges{mcs(imcs)==validMCS}; % SNR points to simulate from MCS

    % 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',10);
        stream.Substream = isnr;
        RandStream.setGlobalStream(stream);

        % Set simulation parameters
        numPacketErrors = 0;
        numPkt = 1; % Index of the transmitted packet

        while numPacketErrors<=maxNumErrors && numPkt<=maxNumPackets
            % Generate a packet waveform
            psdu = randi([0 1],cfgDMG.PSDULength*8,1);
            txWaveform = wlanWaveformGenerator(psdu,cfgDMG);

            % Add delay and trailing zeros
            tx = [zeros(delay,1); txWaveform; zeros(zeroPadding,1)];

            % Transmit through a TGay channel. Reset the channel for a
            % different realization per packet.
            reset(tgayChan);
            chanOut = tgayChan(tx);

            % Add noise
            rx = awgn(chanOut,snr(isnr));

            % Add CFO
            rx = helperFrequencyOffset(rx,fs,freqOffset);

            % Packet detection
            threshold = 0.03; % Good for low SNRs
            pktStartOffset = dmgPacketDetect(rx,0,threshold);
            if isempty(pktStartOffset) % If empty no STF detected; packet error
                numPacketErrors = numPacketErrors+1;
                numPkt = numPkt+1;
                continue; % Go to next loop iteration
            end

            % Frequency offset estimation and correction
            stf = rx(pktStartOffset+(ind.DMGSTF(1):ind.DMGSTF(2)));
            fOffsetEst = dmgCFOEstimate(stf);
            rx = helperFrequencyOffset(rx,fs,-fOffsetEst);

            % Symbol timing and channel estimate
            preamblefield = rx(pktStartOffset+1:pktStartOffset+ind.DMGHeader(2),:);
            [symbolTimingOffset,chanEst] = dmgTimingAndChannelEstimate(preamblefield);
            startOffset = pktStartOffset+symbolTimingOffset;

            % If not enough samples to decode detected data field start,
            % then assume synchronization error and packet error
            if (startOffset+ind.DMGData(2))>size(rx,1)
                numPacketErrors = numPacketErrors+1;
                numPkt = numPkt+1;
                continue; % Go to next loop iteration
            end

            % Noise estimation using the STF as repeating sequence
            stf = rx(pktStartOffset+(ind.DMGSTF(1):ind.DMGSTF(2)));
            nVarEst = dmgSTFNoiseEstimate(stf);

            % Extract data field (ignore first GI)
            rxData = rx(startOffset+((ind.DMGData(1)+Ngi):ind.DMGData(2)));

            % Linear frequency domain equalization
            rxEqDataBlks = dmgSingleCarrierFDE(rxData,chanEst,nVarEst);

            % Unique word phase tracking
            rxEqDataBlks = dmgUniqueWordPhaseTracking(rxEqDataBlks);

            % Discard GI from all blocks
            rxDataSym = rxEqDataBlks(1:end-Ngi,:);

            % Recover the transmitted PSDU from DMG Data field
            dataDecode = wlanDMGDataBitRecover(rxDataSym,nVarEst,cfgDMG);

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

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

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