% Set transmission parameters
numSTAs          = 4; % Number of uplink STAs
STAID            = [1 19 2 4]; % STA association ID, assigned to each associated STA with an AP
multiplexingFlag = [0 1 0 0]; % Signaled in the trigger frame of type NFRP for each STA
feedbackStatus   = [1 0 1 0]; % Resource request information signaled by each STA
startingAID      = 0; % Signaled in the trigger frame of type NFRP

% Set simulation parameters
numPackets = 50; % Number of packets to simulate
snrRange   = -2:2:6; % SNR points (dB)
chanBW     = 'CBW20'; % Channel bandwidth
numTx      = 1; % Assume same number of transmit antennas for all STAs
numRx      = 4; % Number of receive (AP) antennas

% Return the index of the transmitting STAs. Calculate RUToneSetIndex and starting space time stream for all STAs. 
[txSTAIndex,ruToneSetIndexPerSTA,startingSTSNumPerSTA] = heTBNDPMappingParams(chanBW,numSTAs,multiplexingFlag,startingAID,STAID);
numTxSTAs = numel(txSTAIndex); % STAs scheduled to transmit
cfgSTA = cell(1,numTxSTAs);

% Generate the configuration object and set the feedback status property for all STAs
cfgBase = wlanHETBConfig('ChannelBandwidth',chanBW,'NumTransmitAntennas',numTx,'SpatialMapping','Fourier');
cfgNDP = getNDPFeedbackConfiguration(cfgBase);
if numSTAs~=numel(feedbackStatus)
    error('The number of elements in FeedbackStatus must be equal to the number of STAs');
end

for u=1:numTxSTAs
      cfgNDP.RUToneSetIndex = ruToneSetIndexPerSTA(u);
      cfgNDP.StartingSpaceTimeStream = startingSTSNumPerSTA(u);
      cfgNDP.FeedbackStatus = feedbackStatus(txSTAIndex(u));
      cfgSTA{u} = cfgNDP;
end

delayProfile = 'Model-D'; % TGax channel multipath delay profile
tgaxBase = wlanTGaxChannel;
tgaxBase.DelayProfile = delayProfile;
tgaxBase.SampleRate = wlanSampleRate(cfgSTA{1});
tgaxBase.TransmissionDirection = 'Uplink';
tgaxBase.TransmitReceiveDistance = 10;
tgaxBase.ChannelBandwidth = chanBW;
tgaxBase.NumReceiveAntennas = numRx;
tgaxBase.NormalizeChannelOutputs = false;
tgaxBase.PathGainsOutputPort = true;

% A cell array stores the channel objects, one per STA.
tgax = cell(1,numTxSTAs);
for u=1:numTxSTAs
    tgax{u} = clone(tgaxBase);
    tgax{u}.NumTransmitAntennas = numTx;
    tgax{u}.UserIndex = u;
end
chInfo = info(tgaxBase);
chFilterCoefficients = chInfo.ChannelFilterCoefficients; % Channel filter coefficients

% Processing SNR Points
ofdmInfo = wlanHEOFDMInfo('HE-LTF',cfgNDP);
numSNR = numel(snrRange); % Number of SNR points
misdetectionProbability = zeros(numTxSTAs,numSNR);
ind = wlanFieldIndices(cfgNDP); % Same for all STAs

%parfor isnr=1:numSNR % Use 'parfor' to speed up the simulation
for isnr=1:numSNR
    % 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);
    rxFeedbackStatus = zeros(numPackets,numTxSTAs);
    chDelay = zeros(1,numTxSTAs);
      
    for pktIdx=1:numPackets
        rxWaveform = 0;

        % Generate random channel realization for each packet by varying
        % the UserIndex property of the channel. This assumes all STAs
        % have the same number of transmit antennas.
        chPermutations = randperm(numTxSTAs);
        for u=1:numTxSTAs
            % Generate HE TB feedback NDP waveform for each STA
            txSTA = wlanWaveformGenerator([],cfgSTA{u});
            
            % Pass waveform through a random TGax Channel
            channelIdx = chPermutations(u);
            reset(tgax{channelIdx}); % New channel realization
            [rxSTA,h] = tgax{channelIdx}([txSTA; zeros(50,size(txSTA,2))]);
            
            % Perform perfect channel delay estimate to find the start of
            % the packet
            chDelay(u) = channelDelay(h,chFilterCoefficients);

            % Combine uplink waveform from all STAs into one waveform
            rxWaveform = rxWaveform+rxSTA;   
        end

        % Synchronize later in time by using the maximum channel delay
        % between all channels as the start of the packet
        pktOffset = max(chDelay(u)); % Packet start index
        
        % Pass the waveform through AWGN channel. Account for noise
        % energy in unused subcarriers.
        snrVal = snrRange(isnr)-10*log10(ofdmInfo.FFTLength/ofdmInfo.NumTones);
        rxWaveform = awgn(rxWaveform,snrVal);
        
        % Uplink processing (at the AP)
        rxHELTF = rxWaveform(pktOffset+(ind.HELTF(1):ind.HELTF(2)),:);
        heltfDemod = wlanHEDemodulate(rxHELTF,'HE-LTF',chanBW,cfgNDP.GuardInterval,cfgNDP.HELTFType);

        % Recover feedback status for all STAs
        for u=1:numTxSTAs
            rxFeedbackStatus(pktIdx,u) = wlanHETBNDPFeedbackStatus(heltfDemod,cfgSTA{u});
        end
    end
    % Probability of misdetection per STA
    misdetectionProbability(:,isnr) = 1-sum(rxFeedbackStatus==feedbackStatus(txSTAIndex))/numPackets;
    disp(['SNR ' num2str(snrRange(isnr)) ' completed']);
end

markers = 'ox*sd^v><ph+';
numMarkers = numel(markers);
for u=1:numTxSTAs
    semilogy(snrRange,misdetectionProbability(u,:),['-' markers(mod(u-1,numMarkers)+1)]); hold on;
end
xlabel('SNR(dB)');
ylabel('Probability of misdetection');
dataStr = arrayfun(@(x)sprintf('STAID %d',x),STAID(txSTAIndex),'UniformOutput',false);
title(sprintf('%s, Channel %s, %d-by-%d',chanBW,delayProfile,numTx,numRx));
legend(dataStr);
grid on;