NR V2X Sidelink PSCCH and PSSCH Throughput

% Use the PSSCH RMC test configurations as initial parameter templates for the simulation,
% TS 38.101-4 Section 11, V2X requirements and A.6, SL reference measurements channels

% Select PSSCH demodulation test
% All PSSCH RMC below are defined for 20 MHz BW @ 30 kHz SCS (10 PRB subchannels)
% 'R.PSSCH.2-1.1' - 20 PRB (2 subchannels), QPSK,  TCR=308/1024, TDLA30-2700
% 'R.PSSCH.2-1.2' - 20 PRB (2 subchannels), 16QAM, TCR=378/1024, TDLA30-1400
% 'R.PSSCH.2-1.3' - 10 PRB (1 subchannel),  64QAM, TCR=438/1024, TDLA30-180
rmc   = "R.PSSCH.2-1.1";   % RMC and test selection
bwscs = [20 30];   % Release 16 NR sidelink is specified to operate in FR1 at 10, 20, 30, 40 MHz BW

% Get PSCCH/PSSCH channel and sidelink pool parameters
bw = bwscs(1); 
scs = bwscs(2);
[sltransmission,pool,testconfig] = NRSidelinkResourcePool.rmcTestConfiguration(rmc,bw,scs);

% Reconfigure any required pool and channel parameters here for the throughput/BLER simulation
% For example:
% pool.sl_PSFCH_Period_r16 = 0;  % Uncomment to turn off PSFCH symbol reservation in the resource pool (in terms of L3 protocol, no HARQ enabled in the pool)

% Visually crosscheck the PSCCH and PSSCH transmission resources configured for simulation
pool

sltransmission

displayTransmissionResources(pool,sltransmission);

% Define all top-level simulation parameters in a single 'simParameters' structure
simParameters = struct();
simParameters.Pool = pool;
simParameters.SLTransmission = sltransmission;

% From TS 38.101-4 Section 11, V2X requirements tests
% Reconfigure the propagation channel parameters, as required
simParameters.DelayProfile = testconfig.DelayProfile;
simParameters.MaximumDopplerShift = testconfig.MaximumDopplerShift;

% Antenna geometry
simParameters.NTxAnts = sltransmission.NumLayers;
simParameters.NRxAnts = 2;

% Simulation length parameters 
simParameters.SNRIn = testconfig.ReferenceSNR + (-5:5); % Create an SNR range around the test target reference SNR
simParameters.NumSlots = 100;                           % Number of transmission slots to simulate, per SNR point

% HARQ-related parameters
simParameters.NHARQProcesses = max(1,pool.sl_PSFCH_Period_r16);   % Number of independent parallel HARQ processes to simulate (up to 16 processes available (4 bits available for HARQ ID in SCI2))
simParameters.RVSequence = [0];     % Sequence of RV indices for a TB (initial transmission and any potential retransmissions, for full HARQ the sequence [0 2 3 1] is performant)
simParameters.BlindHARQ = 0;        % Enable blind open-loop retransmissions, according to the above specified RV sequence

% General algorithmic parameters
simParameters.PerfectChannelEstimator = 1;  % Use perfect channel estimate; otherwise, use a DM-RS-based estimate
simParameters.NXIDPassthrough = 1;          % Pass the transmission PSSCH scrambling NXID value directly to the receiver

% Simulation diagnostics
simParameters.DisplaySimulationInformation = 0;  % Display information about each transmission event

% Create and configure the propagation channel model object for the simulation
if contains(simParameters.DelayProfile,'CDL','IgnoreCase',true)
    propchannel = nrCDLChannel('DelayProfile',simParameters.DelayProfile);
    % CDL uses the antenna panel array sizing, rather than just the total number of antennas, so conversion is required
    [propchannel.TransmitAntennaArray,simParameters.NTxAnts] = ueArrayGeometry(propchannel.TransmitAntennaArray,simParameters.NTxAnts);
    [propchannel.ReceiveAntennaArray,simParameters.NRxAnts]  = ueArrayGeometry(propchannel.ReceiveAntennaArray,simParameters.NRxAnts);
else
    propchannel = nrTDLChannel('DelayProfile',simParameters.DelayProfile);
    propchannel.NumTransmitAntennas = simParameters.NTxAnts;   % Configure the number of transmit antennas for the transmitting UE
    propchannel.NumReceiveAntennas  = simParameters.NRxAnts;   % Configure the number of receive antennas for the receiving UE
end
propchannel.MaximumDopplerShift = simParameters.MaximumDopplerShift;
waveinfo = getWaveformInfo(simParameters.Pool);            % Get details of the OFDM baseband IQ that will be transmitted in the pool
propchannel.SampleRate = waveinfo.SampleRate;              % Align propagation channel with sample rate of baseband transmission

% Start of main BLER simulation body

% Array to store the maximum SL-SCH throughput for all SNR points
maxThroughput = zeros(length(simParameters.SNRIn),1);
% Array to store the simulation SL-SCH throughput for all SNR points
simThroughput = zeros(length(simParameters.SNRIn),1);

% Record the control/data information block errors for each SNR point
simEtotals = repmat(struct('SCI1',0,'SCI2',0,'SLSCH',0,'TxCount',0),1,length(simParameters.SNRIn));

% Set up SL-SCH transport channel processing
slschtxProto = nrULSCH("MultipleHARQProcesses",true);     % SL-SCH transport channel encoder (same processing as UL-SCH, but with I_BIL = 0 always)
slschrxProto = nrULSCHDecoder("MultipleHARQProcesses",true,'LDPCDecodingAlgorithm','Normalized min-sum');  % SL-SCH decoder

% Use fixed target code rate for all processes and throughout the simulation
slschtxProto.TargetCodeRate = simParameters.SLTransmission.TargetCodeRate;
slschrxProto.TargetCodeRate = simParameters.SLTransmission.TargetCodeRate;

% Loop through each SNR point in overall simulation
for snrIdx = 1:numel(simParameters.SNRIn)        % Serializing the execution of the SNR points
% parfor snrIdx = 1:numel(simParameters.SNRIn)   % Uncomment to parallelize the SNR points using the Parallel Computing Toolbox parfor

    % Accessing properties in objects/structures defined outside of the parfor body
    % may result in broadcast variable warnings, so take a local copy of the entire variable
    simParametersLocal = simParameters;

    % Using a locally created copy of the pool ensures that the pool state (e.g., current slot number)
    % starts the same for each SNR point run
    pool = simParametersLocal.Pool;
    pool.SCSCarrier.NSlot = 0;
    % Take a copy of the PSCCH/PSSCH transmission parameters for ease of reference
    sltransmission = simParametersLocal.SLTransmission;
    
    % Take a copy of the transport encoder/decoder handles to help PCT with classification of variables
    slschtx = slschtxProto;
    slschrx = slschrxProto;
    reset(slschrx);             % Reset SL-SCH decoder at the start of each SNR point (explicitly required for serial 'for' loop)

    % Get the SNR value for this run
    SNRdB = simParametersLocal.SNRIn(snrIdx);

    % Print summary of overall configuration for this run
    fprintf('\nSimulating PSCCH/PSSCH at %gdB SNR (%d/%d) with %s-%d channel (%dx%d) for %d slots \n', ...
                    SNRdB,snrIdx,numel(simParametersLocal.SNRIn), ...
                    simParametersLocal.DelayProfile,simParametersLocal.MaximumDopplerShift,simParametersLocal.NTxAnts,simParametersLocal.NRxAnts, ...
                    simParametersLocal.NumSlots);

    % Reset the random number generator so that each SNR point experiences
    % the same noise realization
    rng(0,'twister');

    % Reset the propagation channel state so that each SNR point experiences
    % the same channel realization
    reset(propchannel);
    
    % Get local copy of info about the baseband waveform to be simulated
    waveinfo = getWaveformInfo(pool);

    % Use a locally created SL-SCH HARQ entity object to ensure independence between SNR points
    % This object controls the use and HARQ scheduling of the SL-SCH transport encoder
    harqEntity = HARQEntitySL(0:simParametersLocal.NHARQProcesses-1,simParametersLocal.RVSequence,1,simParametersLocal.BlindHARQ);          % Up to 16 processes available (4 bits for HARQ ID in SCI2)

    % Record the control/data information block errors for this SNR point
    etotals = simEtotals(snrIdx);
    
    % Simulation loop for a single SNR point
    for nslot=1:simParametersLocal.NumSlots
    
        % Prepare transport/control data for transmission

        % PSFCH slot related adjustments to the PSSCH transmission parameters
        sltransmission.PSFCHOverhead = pool.IsPSFCHSlot;
        if pool.IsPSFCHSlot && ~isempty(sltransmission.NumDMRSPositionsList)
            sltransmission.NumDMRSPositions = min(sltransmission.NumDMRSPositionsList);
        else
            sltransmission.NumDMRSPositions = max(sltransmission.NumDMRSPositionsList);
        end

        % SL-SCH transport data
        if harqEntity.NewData   % Is this a new data opportunity for the current process?
    
            % Calculate the SL-SCH TBS from the transmission and pool parameters
            tbs = getTBS(pool,sltransmission);
    
            % Transmit side SL-SCH preparation
            trBlk = randi([0 1],tbs,1);
            setTransportBlock(slschtx,trBlk,harqEntity.HARQProcessID);
    
            % Receive side SL-SCH preparation
            slschrx.TransportBlockLength = tbs;
            resetSoftBuffer(slschrx,harqEntity.HARQProcessID);          % Reset the process soft buffer if transmitting new data on latest process
        end

        % Update dynamic transmission parameters for the SL-SCH part on PSSCH transmission
        sltransmission.HARQProcessID = harqEntity.HARQProcessID;
        sltransmission.RedundancyVersion = harqEntity.RedundancyVersion;
        sltransmission.NXID = [];   % Ensure there is no external overridding of the PSSCH NXID value in the transmission

        % SCI1 and SCI2 data 
        %
        % 1st-stage SCI information bits for PSCCH
        sci1infobits = randi([0 1],sltransmission.OSCI1,1);
        % 2nd-stage SCI information bits for transmission on PSSCH, multiplexed with SL-SCH data
        sci2infobits = randi([0 1],sltransmission.OSCI2,1);
    
        % Transmit PSCCH/PSSCH in the resource pool (baseband)
        % Input is the separate control information bits (SCI1 and SCI2) and transport channel processing object
        % Output is the slot waveform, plus intermediate processing symbols and resource information
        [basebandiq,sld,slo] = transmitInPool(pool,sltransmission,sci1infobits,sci2infobits,slschtx);  
    
        % Pass baseband waveform through propagation channel
        [basebandiq,pathGains,sampleTimes] = propchannel(basebandiq);
    
        % Add AWGN to the received time domain waveform
        % Normalize noise power by the IFFT size used in OFDM modulation,
        % as the OFDM modulator applies this normalization to the
        % transmitted waveform. Also normalize by the number of receive
        % antennas, as the channel model applies this normalization to the
        % received waveform, by default
        SNR = 10^(SNRdB/10);
        N0 = 1/sqrt(size(basebandiq,2)*double(waveinfo.Nfft)*SNR);
        noise = N0*randn(size(basebandiq),"like",basebandiq);
        basebandiq = basebandiq + noise;
    
        % Parameterize the decoder algorithm
        decalgo = struct('PerfectChannelEstimator',simParametersLocal.PerfectChannelEstimator, ... % Control whether to use perfect or DM-RS-based channel estimation
                         'DecodeSLSCH',~harqEntity.TransmissionSuccess);                           % Control whether to decode the SL-SCH transport block
        if decalgo.PerfectChannelEstimator
            % If perfect channel estimation, then use the value of the path gains
            % provided by the propagation channel
            pathFilters = getPathFilters(propchannel); % Get path filters for perfect channel estimation
            decalgo.estChannelGrid = nrPerfectChannelEstimate(pool.SCSCarrier,pathGains,pathFilters,0,sampleTimes);

            % Get perfect noise estimate (from the noise realization)
            noiseGrid = nrOFDMDemodulate(pool.SCSCarrier,noise);
            decalgo.noiseEst = var(noiseGrid(:));
        end

        % Pass the transmitted PSCCH NXID value to the receiver, if required
        % This operation allows the PSSCH decoding to be independent of the PSCCH decoding success
        if simParametersLocal.NXIDPassthrough
            sltransmission.NXID = sld.NXID;
        end

        % Receive the transmission in the current slot of the pool
        [eset,sldr,slor] = receiveInPool(pool,sltransmission,basebandiq,slschrx,decalgo);
    
        % Accumulate the transmission block error results
        for efn = string(fieldnames(eset))'
            em = eset.(efn);
            if (efn=="SLSCH") && (harqEntity.TransmissionNumber == harqEntity.MinNumTransmissions-1) % Is this the first SL-SCH transmission event for error accounting
                % Record error if this transmission is in error and there have been no passes before on the open process
                % and then scale error count so that it aligns with overall TxCount
                em = (em && ~harqEntity.TransmissionSuccess)*harqEntity.MinNumTransmissions;
            end
            etotals.(efn) = etotals.(efn) + em;  % Accumulate error status into error total
        end
        etotals.TxCount = etotals.TxCount + 1;   % Keep a total event count, needed for the block error rate
        
        % Get the SL-SCH reception result for SL-SCH transmit HARQ management
        blkerr = eset.SLSCH;    
        G = slo.PSSCHG(2);       % G is bit capacity (codeword size) of transmission resources for SL-SCH on PSSCH

        % Update current HARQ process state with the results of the associated SL-SCH/PSSCH transmission
        % and advance to next process
        tbs = length(getTransportBlock(slschtx,harqEntity.HARQProcessID));
        slschrxreport = updateAndAdvance(harqEntity,blkerr,tbs,G);
        
        % Print out transmission status
        if simParameters.DisplaySimulationInformation
            etext = ["passed" "FAILED"];
            fprintf('\n%s %s %s %s', ...
                pad(sprintf('(%3.2f%%)',100*nslot/simParametersLocal.NumSlots),9), ...  % Percentage of simulate complete
                pad(sprintf('NSlot=%d,',pool.SCSCarrier.NSlot),10), ...                 % Absolute slot number
                sprintf('SCI1=%s, SCI2=%s,',etext(1+eset.SCI1),etext(1+eset.SCI2)), ... % SCI error status           
                slschrxreport);                                                         % SL-SCH HARQ information
        end

        % Move on to next available slot in the resource pool
        [pool,nslotsadvanced] = advanceToNextSlot(pool);

    end

    % For the current SNR point, 
    % record the control/data information block errors 
    simEtotals(snrIdx) = etotals;

    % Store values to calculate SL-SCH throughput
    simThroughput(snrIdx) = harqEntity.SuccessfulBitsAllProcesses;
    maxThroughput(snrIdx) = harqEntity.TotalBitsAllProcesses;

end

plotResults(simParameters,simEtotals,simThroughput,maxThroughput);

function [antArray,nants] = ueArrayGeometry(antArray,nUeAnts)
% Turn scalar 'number of antennas' into an antenna panel/panel array dimensions
    % The target panel array is a single panel comprising a row of antenna elements
    % of dual polarization (when the number of antennas specified is greater that 1)
    % 
    % The antenna array is defined in terms of the number of rows and columns of a single antenna panel,
    % the number of polarizations, and the number of rows and columns in the overall array of individual panels
    arraySize = [ceil(nUeAnts/2),1,1+(nUeAnts>1),1,1];
    nants = prod(arraySize);    % Final number of antenna elements in overall array
    antArray.Size = arraySize;
    antArray.Element = 'isotropic';
end

function plotResults(simParameters,simEtotals,simThroughput,maxThroughput)
% Display simulation BLER and throughput results in a new figure
    figure;
    channelestlabel = ["Practical" "Perfect"];
    sgtitle(sprintf('PSCCH/PSSCH Performance (%s, %s-%d, %dx%d)',channelestlabel(1+logical(simParameters.PerfectChannelEstimator)),simParameters.DelayProfile,simParameters.MaximumDopplerShift,simParameters.NTxAnts,simParameters.NRxAnts));
    epn = 1;
    [snrlb, snrub]=bounds(simParameters.SNRIn);
    % Plot the BLER for SCI1, SCI2, and SL-SCH
    for efn = string(fieldnames(simEtotals))'
        subplot(2,2,epn);
        semilogy(simParameters.SNRIn,[simEtotals.(efn)]./[simEtotals.TxCount],'o');
        if snrub > snrlb
            ax = axis; axis([snrlb snrub ax(3:4)]); % Fix the range of the SNR x-axis (semilogy might have rescaled it for SNR points with no errors)
        end
        grid on;
        xlabel('SNR (dB)'); ylabel('Block Error Rate');
        title(strcat(efn,' BLER'));
        epn = epn + 1;
    end
    % Plot throughput for SL-SCH
    plot(simParameters.SNRIn,simThroughput*100./maxThroughput,'o-.')
    if snrub > snrlb
        ax = axis; axis([snrlb snrub ax(3:4)]);     % Fix the range of the SNR x-axis
    end
    grid on;
    xlabel('SNR (dB)'); ylabel('Throughput (%)'); 
    title('SLSCH Throughput');
    % Reannotate the SL-SCH plots with the TCR and modulation
    for i = [3 4]
        subplot(2,2,i)
        subtitle(sprintf('TCR = %.3f, Mod = %s',simParameters.SLTransmission.TargetCodeRate,simParameters.SLTransmission.Modulation));
    end
    % Reannotate the SCI2 plot with the beta offset value
    subplot(2,2,2)
    subtitle(sprintf('TCR = %.3f, Beta = %.3f',simParameters.SLTransmission.TargetCodeRate,simParameters.SLTransmission.BetaSCI2Offset));
end
%#ok<*NBRAK2>