% Create a format configuration object
cfgDMG = wlanDMGConfig;
cfgDMG.MCS = 0; % MCS 0 represents Control PHY
cfgDMG.PSDULength = 256; % PSDULength in bytes

% Define the pulse shaping filter characteristics
pulseShaping = true; % Enable pulse shaping
Nsym = 8;            % Filter span in symbol durations
alpha = 0.5;         % Roll-off factor
osps = 4;            % Output samples per symbol

% Transmit pulse shaping filter
txFilter = comm.RaisedCosineTransmitFilter;
txFilter.RolloffFactor = alpha;
txFilter.FilterSpanInSymbols = Nsym;
txFilter.OutputSamplesPerSymbol = osps;
txFilter.Shape = 'Normal';

% Receive pulse shaping filter
rxFilter = comm.RaisedCosineReceiveFilter;
rxFilter.RolloffFactor = alpha;
rxFilter.DecimationFactor = osps;
rxFilter.InputSamplesPerSymbol = osps;
rxFilter.FilterSpanInSymbols = Nsym;
rxFilter.Shape = 'Normal';

snrVec = -13.5:0.5:-10.5;

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

numSNR = numel(snrVec); % Number of SNR points
packetErrorRate = zeros(numSNR,1);
indices = wlanFieldIndices(cfgDMG);

if ~strcmp(phyType(cfgDMG),'Control')
    error('This example only supports DMG Control PHY simulation');
end

%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.SNR = snrVec(isnr);

    % Noise power
    nVar = 10^(-snrVec(isnr)/10);
    numPacketErrors = 0;
    numPkt = 1; % Index of packet transmitted

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

        % Transmitter filtering
        if pulseShaping
            % Append zero to compensate for filter group delay
            tx = txFilter([tx; zeros(Nsym,1)]);
            reset(txFilter);
        end

        % Add noise
        rx = awgnChannel(tx);

        % Receiver filtering
        if pulseShaping
            rx = rxFilter(rx);
            reset(rxFilter);
        end

        % Synchronize
        % The received signal is synchronized to the start of the packet by
        % compensating for a known delay due the spectral shaping filters
        if pulseShaping
            offset = Nsym;
        else
            offset = 0; %#ok<UNRCH>
        end

        % Process header and data field together
        rxHeaderDataField = rx(offset+(indices.DMGHeader(1):indices.DMGData(2)));

        % Apply pi/2 de-rotation and de-spread the received signal
        [rxSym,SF] = dmgControlDespread(rxHeaderDataField);

        % Recover the transmitted PSDU from DMG Data field. Scale the noise
        % power by the spreading factor
        dataDecode = wlanDMGDataBitRecover(rxSym,nVar/SF,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(isnr) = numPacketErrors/(numPkt-1);
    disp(['SNR ' num2str(snrVec(isnr))...
          ' completed after ' num2str(numPkt-1) ' packets,'...
          ' PER: ' num2str(packetErrorRate(isnr))]);
end

figure;
semilogy(snrVec,packetErrorRate,'-o');
grid on;
xlabel('SNR (dB)');
ylabel('PER');
legend('MCS 0');
title('PER for DMG Control-PHY, AWGN channel');

function [y,SF] = dmgControlDespread(rx)
    SF = 32; % Spreading factor
    dataField = rx.*exp(-1i*pi/2*(0:size(rx,1)-1).'); % De-rotate symbols
    Ga = wlanGolaySequence(SF); % Generate Golay sequence
    y = (reshape(dataField,SF,length(dataField)/SF)'*Ga)/SF;
end