cfgVHT = wlanVHTConfig;             % Create VHT transmission configuration
cfgVHT.ChannelBandwidth = 'CBW160'; % Bandwidth
cfgVHT.MCS = 7;                     % 64-QAM, rate 5/6
NumReceiveAntennas = 2;             % Number of receive antennas

cfgVHT.APEPLength = 4096; % Bytes
cfgVHT.STBC = false;
cfgVHT.NumTransmitAntennas = NumReceiveAntennas;
cfgVHT.NumSpaceTimeStreams = NumReceiveAntennas;
cfgVHT.SpatialMapping = 'Direct';
cfgVHT.GuardInterval = 'Long';

testInputLevelOffsets = [-10 -9 -8 -7]; % dB

maxNumErrors = 20;
maxNumPackets = 200;

% Receiver minimum input level sensitivity for 20 MHz, Table 22-25. The
% sensitivity increases by 3dB for double the bandwidth.
rxMinSensitivityTable = [-82 -79 -77 -74 -70 -66 -65 -64 -59 -57]; % dBm

% Get minimum input sensitivity given MCS and bandwidth
fs = wlanSampleRate(cfgVHT);  % Baseband sampling rate (Hz)
B = floor(10*log10((fs/20e6))); % Scalar for bandwidth
rxMinSensitivity = rxMinSensitivityTable(cfgVHT.MCS+1)+B; % dBm
disp(['Minimum sensitivity for MCS' num2str(cfgVHT.MCS) ', ' ...
    num2str(fs/1e6) ' MHz: ' num2str(rxMinSensitivity,'%2.1f') ' dBm'])

testInputLevels = rxMinSensitivity+testInputLevelOffsets; % dBm

A = 10.^((testInputLevels-30)/20);      % Voltage gain (attenuation)
A = A*sqrt(cfgVHT.NumTransmitAntennas); % Account for generator scaling

NF = 6;         % Noise figure (dB)
T = 290;        % Ambient temperature (K)
BW = fs;        % Bandwidth (Hz)
k = 1.3806e-23; % Boltzmann constant (J/K)
noiseFloor = 10*log10(k*T*BW)+NF; % dB
disp(['Receiver noise floor: ' num2str(noiseFloor+30,'%2.1f') ' dBm'])

awgnChannel = comm.AWGNChannel('NoiseMethod','Variance', ...
    'Variance',10^(noiseFloor/10));

ind = wlanFieldIndices(cfgVHT); % For accessing fields within the packet
chanBW = cfgVHT.ChannelBandwidth;
rng(0);  % Set random state for repeatability

agc = comm.AGC; % Automatic gain control

S = numel(testInputLevels);
packetErrorRate = zeros(S,1);
rxAntennaPower = zeros(S,1);
for i=1:S
    disp(['Simulating ' num2str(testInputLevels(i),'%2.1f') ...
        ' dBm input level...']);

    % Loop to simulate multiple packets
    numPacketErrors = 0;
    measuredPower = zeros(maxNumPackets,1); % Average power per antenna
    numPkt = 1; % Index of packet transmitted
    while numPacketErrors<=maxNumErrors && numPkt<=maxNumPackets
        % Generate a packet waveform
        txPSDU = randi([0 1],cfgVHT.PSDULength*8,1); % PSDULength in bytes
        tx = wlanWaveformGenerator(txPSDU,cfgVHT);

        % Scale input signal to desired level
        rx = tx.*A(i);

        % Measure the average power at the antenna connector in Watts
        measuredPower(numPkt) = mean(mean(rx.*conj(rx)));

        % Add noise floor at receiver
        rx = awgnChannel(rx);

        % Pass each channel through AGC
        for ic = 1:size(rx,2)
            rx(:,ic) = agc(rx(:,ic));
            reset(agc);
        end

        % 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 out of a reasonable range (>50 samples);
        % 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
        lltf = rx(pktOffset+(ind.LLTF(1):ind.LLTF(2)),:);
        fineFreqOff = wlanFineCFOEstimate(lltf,chanBW);
        rx = helperFrequencyOffset(rx,fs,-fineFreqOff);

        % Extract VHT-LTF samples from the waveform, demodulate and perform
        % channel estimation
        vhtltf = rx(pktOffset+(ind.VHTLTF(1):ind.VHTLTF(2)),:);
        vhtltfDemod = wlanVHTLTFDemodulate(vhtltf,cfgVHT);
        chanEst = wlanVHTLTFChannelEstimate(vhtltfDemod,cfgVHT);

        % Get single stream channel estimate
        chanEstSSPilots = vhtSingleStreamChannelEstimate(vhtltfDemod,cfgVHT);

        % Extract VHT Data samples from the waveform
        vhtdata = rx(pktOffset+(ind.VHTData(1):ind.VHTData(2)),:);

        % Estimate the noise power in VHT data field
        nEstVHT = vhtNoiseEstimate(vhtdata,chanEstSSPilots,cfgVHT);

        % Recover the transmitted PSDU in VHT Data
        rxPSDU = wlanVHTDataRecover(vhtdata,chanEst,nEstVHT,cfgVHT);

        % 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 input level point
    packetErrorRate(i) = numPacketErrors/(numPkt-1);
    disp(['  Completed after ' ...
        num2str(numPkt-1) ' packets, PER: ' ...
        num2str(packetErrorRate(i))]);

    % Calculate average input power per antenna
    rxAntennaPower(i) = 10*log10(mean(measuredPower(1:(numPkt-1))))+30;
    disp(['  Measured antenna connector power: ' ...
        num2str(rxAntennaPower(i),'%2.1f') ' dBm']);
end

figure
semilogy(rxAntennaPower,packetErrorRate,'o-')
hold on
semilogy(rxMinSensitivity,0.1,'rx')
currentxlim = xlim(gca);
xlim([currentxlim(1) currentxlim(2)+1])
grid on
xlabel('Measured power per antenna connector (dBm)');
ylabel('PER');
legend('Simulated PER performance','Maximum PER at minimum sensitivity');
title(sprintf(['Minimum Input Sensitivity Test: MCS%d, %d MHz, ' ...
    '%d Antennas'],cfgVHT.MCS,fs/1e6,cfgVHT.NumTransmitAntennas))