calibrate = true;             % To execute Part A calibration test
systemLevelSimulation = true; % To execute Part B system-level simulation
showScenarioPlot = true;      % To show dynamic simulation plotting updates

PHYParameters = struct;
PHYParameters.TxPower = 20;      % Transmitter power in dBm
PHYParameters.TxGain = 0;        % Transmitter antenna gain in dBi
PHYParameters.RxGain = -2;       % Receiver antenna gain in dBi
PHYParameters.NoiseFigure = 7;   % Receiver noise figure in dB
PHYParameters.NumTxAntennas = 1; % Number of transmitter antennas
PHYParameters.NumRxAntennas = 1; % Number of receiver antennas
PHYParameters.ChannelBandwidth = 'CBW80'; % Bandwidth of system
PHYParameters.TransmitterFrequency = 5e9; % Transmitter frequency in Hz

MACParameters = struct;
MACParameters.NumChannels = 3; % Number of non-overlapping channels
MACParameters.CCALevel = -70;  % Transmission threshold in CCA algorithm (dBm)

% Number of Rooms in [x,y,z] directions
ScenarioParameters = struct;
ScenarioParameters.BuildingLayout = [10 2 5];

% Size of each room in meters [x,y,z]
ScenarioParameters.RoomSize = [10 10 3];

% Number of receivers per room. Note that only one receiver (STA) can be
% active at any given time.
ScenarioParameters.NumRxPerRoom = 1;

SimParameters = struct;
SimParameters.NumDrops = 3;
SimParameters.NumTxEventsPerDrop = 2;

% Total number of transmitters, assuming one transmitter (tx) per room
numTx = prod(ScenarioParameters.BuildingLayout);

% Create transmitter sites with and isotropic antenna element
roomNames = strings(1,numTx);
for siteInd = 1:numTx
    roomNames(siteInd) = "Room " + siteInd;
end
txs = txsite('cartesian','Name',roomNames,...
    'TransmitterFrequency',PHYParameters.TransmitterFrequency, ...
    'TransmitterPower',10.^((PHYParameters.TxPower+PHYParameters.TxGain-30)/10),...
    'Antenna','isotropic');

% Total number of receivers, assuming one transmitter (tx) per room
numRx = numTx*ScenarioParameters.NumRxPerRoom;

% Create receiver sites
roomNames = strings(1,numRx);
for siteInd = 1:numRx
    roomNames(siteInd) = "Room " + (mod(siteInd-1,numTx)+1) + "-" + ceil(siteInd/numTx);
end
rxs = rxsite('cartesian','Name',roomNames,'Antenna','isotropic');

% Receiver noise power in dBm
T = 290; % Temperature (Kelvin)
k = physconst('Boltzmann'); % Boltzmann constant
% Sample rate (Hz)
fs = wlanSampleRate(wlanHESUConfig('ChannelBandwidth',PHYParameters.ChannelBandwidth));
rxNoisePower = 10*log10(k*T*fs)+30+PHYParameters.NoiseFigure;

seed = rng(6); % Seed random number generator and store state

if showScenarioPlot
    hGrid = tgaxBuildResidentialGrid(ScenarioParameters.RoomSize,ScenarioParameters.BuildingLayout, ...
        numTx,numRx,MACParameters.NumChannels);
end

if calibrate

fprintf('Running calibration ...\n');

% Pre-allocate output
output = struct;
output.sinr = zeros(SimParameters.NumDrops,numTx); % for storing SINR values
for drop = 1:SimParameters.NumDrops
    % Drop receivers in each room
    [association,txChannels,rxChannels,txPositions,rxPositions] = tgaxDropNodes( ...
        txs,rxs,ScenarioParameters,MACParameters.NumChannels);

    % All transmitters active
    activeTx = true(numTx,1);

    % Only pick one receiver per Room
    rxAlloc = randi([1 ScenarioParameters.NumRxPerRoom],numTx,1);
    activeRx = reshape(rxAlloc==1:ScenarioParameters.NumRxPerRoom,[],1);

    % Generate propagation model
    propModel = TGaxResidential('roomSize',ScenarioParameters.RoomSize);

    % Get the index of the transmitter for each receiver
    tnum = repmat((1:numTx),1,numRx/numTx);

    % SINR calculation - loop over each non-overlapping channel
    numChannels = numel(unique(txChannels));
    for k = 1:numChannels
        % Use kth non-overlapping channel
        tind = txChannels == k;
        rind = false(size(activeRx));
        rind(activeRx) = rxChannels(activeRx) == k;
        % Get the index of the transmitter of interest for each active receiver
        tsigind = tnum(rind);

        % Calculate SNR
        output.sinr(drop,tind) = sinr(rxs(rind),txs(tind),...
            'ReceiverGain',PHYParameters.RxGain,...
            'ReceiverNoisePower',rxNoisePower,...
            'PropagationModel',propModel,...
            'SignalSource',txs(tsigind));
    end

    % Plot nodes and links
    if showScenarioPlot % update plot data
        mask = txChannels==rxChannels';
        tgaxUpdatePlot(hGrid,txPositions,rxPositions,activeTx,activeRx,mask,txChannels,rxChannels, ...
            sprintf('Box 1 Test 2 "downlink only" calibration, drop #%d/%d',drop,SimParameters.NumDrops));
    end
end

% Plot the CDF of SINR and compare with calibration curves
tgaxCalibrationCDF(output.sinr,'SS1Box1Test2','Long-term Radio Characteristics');

fprintf('Calibration complete \n')

end

if systemLevelSimulation

% Pre-allocate outputs
output.numPkts = zeros(numRx,1);
output.numPktErrors = zeros(numRx,1);
output.sinrMeas = nan(numRx,SimParameters.NumTxEventsPerDrop,SimParameters.NumDrops);
output.sinrEst = nan(numRx,SimParameters.NumTxEventsPerDrop,SimParameters.NumDrops);
output.pktErrorRate = 0;

% For each possible transmitter create a waveform configuration. In this
% example the link and radio parameters are the same for all nodes.
cfgHEBase = wlanHESUConfig;
cfgHEBase.ChannelBandwidth = PHYParameters.ChannelBandwidth; % Channel bandwidth
cfgHEBase.NumTransmitAntennas = PHYParameters.NumTxAntennas; % Number of transmit antennas
cfgHEBase.SpatialMapping = 'Fourier';  % Spatial mapping matrix
cfgHEBase.NumSpaceTimeStreams = 1;     % Number of space-time streams
cfgHEBase.GuardInterval = 0.8;         % Guard interval duration
cfgHEBase.HELTFType = 4;               % HE-LTF compression mode
cfgHEBase.APEPLength = 1e3;            % Payload length in bytes
cfgHEBase.ChannelCoding = 'LDPC';      % Channel coding
cfgHEBase.MCS = 3;                     % Modulation and coding scheme
cfgHE = cell(numTx,1);
for txidx = 1:numTx
    cfgHE{txidx} = cfgHEBase;
end

fprintf('Running System Level Simulation ...\n')
for drop = 1:SimParameters.NumDrops
    fprintf('  Running drop #%d/%d ...\n',drop,SimParameters.NumDrops);

    % Drop receivers in each room
    [association,txChannels,rxChannels,txPositions,rxPositions] = tgaxDropNodes( ...
        txs,rxs,ScenarioParameters,MACParameters.NumChannels);

    % Generate propagation model
    propModel = TGaxResidential('roomSize',ScenarioParameters.RoomSize);

    % Calculate signal strength for all links
    signalStrength = sigstrength(rxs,txs,propModel,'Type','power',...
        'ReceiverGain',PHYParameters.RxGain); % all signal strengths in dBm

    % Threshold signals below noise level to reduce simulation time
    signalStrength(signalStrength < rxNoisePower) = -Inf;

    % Threshold signals that are not on same non-overlapping channel
    signalStrength(~(txChannels == rxChannels')) = -Inf;

    % Mask the transmitter-receiver links that are non-negligible to
    % simulate and get the linear indices
    nonnegligibleMask = signalStrength > -Inf;

    % Reset the non-negligible channels to create a new realization for the
    % current drop
    nonnegligibleIdx = find(nonnegligibleMask)';

    % For each possible active link in a drop create a channel
    % configuration. In this example the link and radio parameters are the
    % same for all nodes.
    tgaxChan = cell(numel(nonnegligibleIdx),1);
    for i = 1:numel(nonnegligibleIdx)
        % Index of transmitter for a given link
        txIdx = mod(nonnegligibleIdx(i)-1,numTx)+1;
        % Channel configuration. The channel realization for each link is
        % different as the global random stream is used.
        tgaxChanBase = wlanTGaxChannel;
        tgaxChanBase.DelayProfile = 'Model-D';
        tgaxChanBase.NumTransmitAntennas = cfgHE{txIdx}.NumTransmitAntennas;
        tgaxChanBase.NumReceiveAntennas = PHYParameters.NumRxAntennas;
        tgaxChanBase.TransmitReceiveDistance = 10; % Distance in meters for NLOS
        tgaxChanBase.ChannelBandwidth = cfgHE{txIdx}.ChannelBandwidth;
        tgaxChanBase.LargeScaleFadingEffect = 'None';
        tgaxChanBase.EnvironmentalSpeed = 0; % m/s, stationary
        tgaxChanBase.SampleRate = fs;
        tgaxChanBase.NormalizeChannelOutputs = false;

        % Store in cell array and reset the channel to generate a new
        % response
        tgaxChan{i} = tgaxChanBase;
        reset(tgaxChan{i});
    end

    for txevent = 1:SimParameters.NumTxEventsPerDrop
        fprintf('    Running transmission event #%d/%d ...\n',txevent,SimParameters.NumTxEventsPerDrop);

        % Determine active transmitters and receivers with Clear Channel Assessment
        [activeTx,activeRx] = tgaxCCA(signalStrength,MACParameters.CCALevel);

        % Plot scenario and links
        if showScenarioPlot
            tgaxUpdatePlot(hGrid,txPositions,rxPositions,activeTx,activeRx,nonnegligibleMask,txChannels,rxChannels, ...
                sprintf('PHY System-Level Simulation, Drop #%d/%d, Transmission Event #%d/%d', ...
                drop,SimParameters.NumDrops,txevent,SimParameters.NumTxEventsPerDrop));
        end

        % Extract elements for active links using activeTx and activeRx
        cfgHEActive = cfgHE(activeTx);
        associationActive = association(activeTx,activeRx);
        nonnegligibleMaskActive = nonnegligibleMask(activeTx,activeRx);
        signalStrengthActive = signalStrength(activeTx,activeRx);

        % Create array containing active channels
        tgaxChanActive = cell(size(associationActive));
        matchIdx = nonnegligibleIdx==find(activeTx&activeRx');
        tgaxChanActive(nonnegligibleMaskActive) = tgaxChan(any(matchIdx,1));

        % Generate a waveform for each non-negligible active link and
        % combine waveforms for each receiver
        [rxWavs,txPSDUActive,signalPower,interfPower] = tgaxGenerateRxWaveforms( ...
            cfgHEActive,tgaxChanActive,nonnegligibleMaskActive,signalStrengthActive,associationActive);

        % Run PHY link simulation for each link and determine if the packet
        % has been decoded successfully. The estimated interference power
        % is passed to the receiver in the place of an interference power
        % measurement algorithm.
        numActiveRxs = sum(activeRx);
        pktError = false(numActiveRxs,1);
        sinrMeas = nan(numActiveRxs,1);
        for rxIdx = 1:numActiveRxs
            [pktError(rxIdx),sinrMeas(rxIdx)] = tgaxModelPHYLink( ...
                rxWavs{rxIdx},cfgHEActive{rxIdx},rxNoisePower,interfPower(rxIdx),txPSDUActive{rxIdx});
        end

        % Store output for active receivers
        output.numPktErrors(activeRx) = output.numPktErrors(activeRx)+pktError;
        output.numPkts(activeRx) = output.numPkts(activeRx) + 1;
        output.sinrMeas(activeRx,txevent,drop) = sinrMeas;

        % Calculate the expected SINR at each receiver
        sinrEst = 10*log10(signalPower./(interfPower+10^((rxNoisePower-30)/10)));
        output.sinrEst(activeRx,txevent,drop) = sinrEst;
    end
end

% Calculate average packet error rate
output.pktErrorRate = sum(output.numPktErrors)/sum(output.numPkts);

disp('Simulation complete')
disp(['Average packet error rate for transmitters: ' num2str(output.pktErrorRate)]);

end

rng(seed); % Restore random state