6G Link-Level Simulation

simParameters = struct();       % Simulation parameters structure
simParameters.NFrames = 2;      % Number of 10 ms frames
simParameters.SNRdB = -10:2:-6; % SNR range (dB)

simParameters.enableParallelism = true;

simParameters.PerfectChannelEstimator = false;

simParameters.DisplaySimulationInformation = true;

simParameters.DisplayDiagnostics = false;

% Set carrier parameters
simParameters.Carrier = pre6GCarrierConfig;                       % Carrier resource grid configuration
simParameters.Carrier.NSizeGrid = 330;                            % Bandwidth in number of resource blocks
simParameters.Carrier.SubcarrierSpacing = 120;                    % Subcarrier spacing

% Set PDSCH parameters
simParameters.PDSCH = pre6GPDSCHConfig;                           % PDSCH definition for all PDSCH transmissions in the BLER simulation

% Define PDSCH time-frequency resource allocation per slot to be full grid (single full grid BWP) and number of layers
simParameters.PDSCH.PRBSet = 0:simParameters.Carrier.NSizeGrid-1;                % PDSCH PRB allocation
simParameters.PDSCH.SymbolAllocation = [0,simParameters.Carrier.SymbolsPerSlot]; % Starting symbol and number of symbols of each PDSCH allocation
simParameters.PDSCH.NumLayers = 1;                                               % Number of PDSCH transmission layers

% This structure is to hold additional simulation parameters for the DL-SCH and PDSCH
simParameters.PDSCHExtension = struct();             

% Define codeword modulation and target coding rate
% The number of codewords is directly dependent on the number of layers so ensure that layers are set first before getting the codeword number
if simParameters.PDSCH.NumCodewords > 1                           % Multicodeword transmission (when number of layers is > 4)
    simParameters.PDSCH.Modulation = {'16QAM','16QAM'};           % 'QPSK', '16QAM', '64QAM', '256QAM', '1024QAM'
    simParameters.PDSCHExtension.TargetCodeRate = [490 490]/1024; % Code rate used to calculate transport block sizes
else
    simParameters.PDSCH.Modulation = '16QAM';                     % 'QPSK', '16QAM', '64QAM', '256QAM', '1024QAM', '4096QAM'
    simParameters.PDSCHExtension.TargetCodeRate = 490/1024;       % Code rate used to calculate transport block sizes
end

% Disable PT-RS
simParameters.PDSCH.EnablePTRS = false;

% PDSCH PRB bundling (TS 38.214 Section 5.1.2.3)
simParameters.PDSCHExtension.PRGBundleSize = [];                  % Any positive power of 2, or [] to signify "wideband"

% HARQ process parameters
simParameters.PDSCHExtension.NHARQProcesses = 16;                 % Number of parallel HARQ processes to use
simParameters.PDSCHExtension.EnableHARQ = true;                   % Enable retransmissions for each process, using RV sequence [0,2,3,1]

% LDPC decoder parameters
simParameters.PDSCHExtension.LDPCDecodingAlgorithm = 'Normalized min-sum';
simParameters.PDSCHExtension.MaximumLDPCIterationCount = 20;

% Number of antennas
simParameters.NTxAnts = 32;                                       % Number of antennas (1,2,4,8,16,32,64,128,256,512,1024) >= NumLayers
simParameters.NRxAnts = 2;

% Define the general CDL propagation channel parameters
simParameters.DelayProfile = 'CDL-A';   
simParameters.DelaySpread = 10e-9;
simParameters.MaximumDopplerShift = 70;

% Cross-check the PDSCH configuration parameters against the channel geometry 
validateParameters(simParameters);

if (simParameters.enableParallelism && canUseParallelPool)
    pool = gcp; % create parallel pool, requires PCT
    numWorkers = pool.NumWorkers;
    maxNumWorkers = pool.NumWorkers;
else
    if (~canUseParallelPool && simParameters.enableParallelism)
        warning("Ignoring the value of enableParallelism ("+simParameters.enableParallelism+")"+newline+ ...
            "The simulation will run using serial execution."+newline+"You need a license of Parallel Computing Toolbox to use parallelism.")
    end
    numWorkers = 1;    % No parallelism
    maxNumWorkers = 0; % Used to convert the parfor-loop into a for-loop
end

str1 = RandStream('Threefry','Seed',1);
constantStream = parallel.pool.Constant(str1);

numSlotsPerWorker = ceil((simParameters.NFrames*simParameters.Carrier.SlotsPerFrame)/numWorkers);
disp("Parallelism: "+simParameters.enableParallelism)

disp("Number of workers: "+numWorkers)

disp("Number of slots per worker: "+numSlotsPerWorker)

disp("Total number of frames: "+(numSlotsPerWorker*numWorkers)/simParameters.Carrier.SlotsPerFrame)

% Results storage
result = struct(NumSlots=0,NumBits=0,NumCorrectBits=0);
results = repmat(result,numWorkers,numel(simParameters.SNRdB));

% Parallel processing, worker parfor-loop
parfor (pforIdx = 1:numWorkers,maxNumWorkers)     
    % Set random streams to ensure repeatability
    % Use substreams in the generator so each worker uses mutually independent streams
    stream = constantStream.Value;      % Extract the stream from the Constant
    stream.Substream = pforIdx;         % Set substream value = parfor index
    RandStream.setGlobalStream(stream); % Set global stream per worker

    % Per worker processing
    results(pforIdx,:) = pdschLink(simParameters,numSlotsPerWorker,pforIdx);
end

[throughput,throughputMbps,summaryTable] = processResults(simParameters,results);
disp(summaryTable)

figure;
plot(simParameters.SNRdB,throughput,'o-.')
xlabel('SNR (dB)'); ylabel('Throughput (%)'); grid on;
title(sprintf('%s (%dx%d) / NRB=%d / SCS=%dkHz', ...
              simParameters.DelayProfile,simParameters.NTxAnts,simParameters.NRxAnts, ...
              simParameters.Carrier.NSizeGrid,simParameters.Carrier.SubcarrierSpacing));

function [throughput,throughputMbps,summaryTable] = processResults(simParameters,results)
% Process multi-worker and multi-SNR results

    numSNRPts = size(results,2);
    
    totalSimulatedSlots = sum(reshape([results(:).NumSlots].',[],numSNRPts),1);
    totalSimulatedBits = sum(reshape([results(:).NumBits].',[],numSNRPts),1);
    totalCorrectBits = sum(reshape([results(:).NumCorrectBits].',[],numSNRPts),1);
    totalSimulatedFrames = totalSimulatedSlots/simParameters.Carrier.SlotsPerFrame;
    
    % Throughput results calculation
    throughput = 100*(totalCorrectBits./totalSimulatedBits);
    throughputMbps = 1e-6*totalCorrectBits/(simParameters.NFrames*10e-3);
    summaryTable = table(simParameters.SNRdB.',totalSimulatedBits.',totalSimulatedSlots.', ...
        totalSimulatedFrames.',throughput.',throughputMbps.');
    summaryTable.Properties.VariableNames = ["SNR" "Simulated bits" "Number of Tr Blocks" ...
        "Number of frames" "Throughput (%)" "Throughput (Mbps)"];

end

function resultsPerWorker = pdschLink(simParameters,totalNumSlots,workerId)
% PDSCH link simulation
    
    % Take copies of channel-level parameters to simplify subsequent parameter referencing 
    carrier = simParameters.Carrier;
    pdsch = simParameters.PDSCH;
    pdschextra = simParameters.PDSCHExtension;

    % Results storage
    result = struct(NumSlots=0,NumBits=0,NumCorrectBits=0);
    resultsPerWorker = repmat(result,1,numel(simParameters.SNRdB));

    % Create DL-SCH encoder/decoder
    [encodeDLSCH,decodeDLSCH] = dlschEncoderDecoder(pdschextra);

    % OFDM waveform information
    ofdmInfo = hpre6GOFDMInfo(carrier);
        
    % Create CDL channel
    channel = nrCDLChannel;
    channel = hArrayGeometry(channel,simParameters.NTxAnts,simParameters.NRxAnts);
    nTxAnts = prod(channel.TransmitAntennaArray.Size);
    nRxAnts = prod(channel.ReceiveAntennaArray.Size);
    channel.DelayProfile = simParameters.DelayProfile;
    channel.DelaySpread = simParameters.DelaySpread;
    channel.MaximumDopplerShift = simParameters.MaximumDopplerShift;
    channel.SampleRate = ofdmInfo.SampleRate;
    % New seed for each worker, but the same for each SNR point so they all
    % experience the same channel realization.
    channel.Seed = randi([0 2^32-1]);

    chInfo = info(channel);
    maxChDelay = chInfo.MaximumChannelDelay;

    % Set up redundancy version (RV) sequence for all HARQ processes
    if simParameters.PDSCHExtension.EnableHARQ        
        rvSeq = [0 2 3 1];
    else
        % HARQ disabled - single transmission with RV=0, no retransmissions
        rvSeq = 0;
    end

    % for all SNR points
    for snrIdx = 1:length(simParameters.SNRdB)

        % Noise power calculation
        SNR = 10^(simParameters.SNRdB(snrIdx)/10); % Calculate linear noise gain
        N0 = 1/sqrt(double(ofdmInfo.Nfft)*SNR*nRxAnts);
        % Get noise power per resource element (RE) from noise power in the
        % time domain (N0^2)
        nPowerPerRE = N0^2*ofdmInfo.Nfft;

        % Reset the channel and DL-SCH decoder at the start of each SNR simulation
        reset(channel);
        reset(decodeDLSCH);        

        % Specify the fixed order in which we cycle through the HARQ process IDs
        harqSequence = 0:pdschextra.NHARQProcesses-1;

        % Initialize the state of all HARQ processes
        harqEntity = HARQEntity(harqSequence,rvSeq,pdsch.NumCodewords);

        % Obtain a precoding matrix (wtx) to be used in the transmission of the
        % first transport block
        estChannelGrid = getInitialChannelEstimate(carrier,nTxAnts,channel);    
        wtx = hSVDPrecoders(carrier,pdsch,estChannelGrid,pdschextra.PRGBundleSize);

        %  Progress when parallel processing is enabled
        if (simParameters.enableParallelism && workerId==1)
            fprintf('Simulating SNR=%.2f dB, progress: \n0%% \n',simParameters.SNRdB(snrIdx))
        end

        % Process all the slots per worker
        for nSlot = 0:totalNumSlots-1

            % New slot number
            carrier.NSlot = nSlot;

            % Calculate the transport block sizes for the transmission in the slot
            [pdschIndices,pdschIndicesInfo] = hpre6GPDSCHIndices(carrier,pdsch);
            trBlkSizes = getTBS(pdsch.Modulation,pdsch.NumLayers,numel(pdsch.PRBSet),pdschIndicesInfo.NREPerPRB,pdschextra.TargetCodeRate);

            % Generate new data and DL-SCH encode
            codedTrBlocks = getDLSCHCodeword(encodeDLSCH,trBlkSizes,pdsch.Modulation,pdsch.NumLayers,pdschIndicesInfo.G,harqEntity);
        
            % PDSCH modulation of codeword(s), MIMO precoding and OFDM
            [txWaveform,pdschSymbols]= hPDSCHTransmit(carrier,pdsch,codedTrBlocks,wtx);

            % Pass data through channel model
            txWaveform = [txWaveform; zeros(maxChDelay,size(txWaveform,2))];
            [rxWaveform,pathGains,sampleTimes] = channel(txWaveform);

            % Add noise
            noise = N0*randn(size(rxWaveform),"like",rxWaveform);
            rxWaveform = rxWaveform + noise;

            % Synchronization, OFDM demodulation, channel estimation,
            % equalization, and PDSCH demodulation            
            pathFilters = getPathFilters(channel);
            
            perfEstConfig = perfectEstimatorConfig(pathGains,sampleTimes,pathFilters,nPowerPerRE,simParameters.PerfectChannelEstimator);
            [dlschLLRs,wtx,pdschEq] = hPDSCHReceive(carrier,pdsch,pdschextra,rxWaveform,wtx,perfEstConfig);

            % Display EVM per layer, per slot and per RB
            if (simParameters.DisplayDiagnostics)
                gridSize = [carrier.NSizeGrid*12 carrier.SymbolsPerSlot nTxAnts];
                plotLayerEVM(totalNumSlots,nSlot,pdsch,gridSize,pdschIndices,pdschSymbols,pdschEq,simParameters.SNRdB(snrIdx));
            end

            % Decode the DL-SCH transport channel

            % If new data because of previous RV sequence time out then flush decoder soft buffer explicitly
            for cwIdx = 1:pdsch.NumCodewords
                if harqEntity.NewData(cwIdx) && harqEntity.SequenceTimeout(cwIdx)
                    resetSoftBuffer(decodeDLSCH,cwIdx-1,harqEntity.HARQProcessID);
                end
            end
            decodeDLSCH.TransportBlockLength = trBlkSizes;
            blkerr = getTransportBlockCRC(decodeDLSCH,dlschLLRs,pdsch,harqEntity); 

            % Update current process with CRC error and advance to next process
            procstatus = updateAndAdvance(harqEntity,blkerr,trBlkSizes,pdschIndicesInfo.G);
            if (simParameters.DisplaySimulationInformation && ~simParameters.enableParallelism)
                fprintf('(%3.2f%%), SNR=%.2f dB, NSlot=%d, %s\n',100*(nSlot+1)/totalNumSlots,simParameters.SNRdB(snrIdx),nSlot,procstatus);
            elseif (simParameters.enableParallelism && workerId==1 && ~mod((nSlot+1),ceil((10*totalNumSlots)/100))) % Progress when parallel processing is enabled
                % Update progress with approximately 10% steps
                fprintf('%3.2f%% \n',100*(nSlot+1)/totalNumSlots)
            end

            % SNR point simulation results            
            resultsPerWorker(snrIdx).NumSlots = resultsPerWorker(snrIdx).NumSlots+1;
            resultsPerWorker(snrIdx).NumBits = resultsPerWorker(snrIdx).NumBits+sum(trBlkSizes);
            resultsPerWorker(snrIdx).NumCorrectBits = resultsPerWorker(snrIdx).NumCorrectBits+sum(~blkerr .* trBlkSizes);

        end % for nSlot = 0:totalNumSlots

    end % for all SNR points
end

function [encodeDLSCH,decodeDLSCH] = dlschEncoderDecoder(PDSCHExtension)
% Create and parameterize the DL-SCH encoder and decoder objects
    
    % Create DL-SCH encoder object
    encodeDLSCH = nrDLSCH;
    encodeDLSCH.MultipleHARQProcesses = true;
    encodeDLSCH.TargetCodeRate = PDSCHExtension.TargetCodeRate;
    
    % Create DL-SCH decoder object
    decodeDLSCH = nrDLSCHDecoder;
    decodeDLSCH.MultipleHARQProcesses = true;
    decodeDLSCH.TargetCodeRate = PDSCHExtension.TargetCodeRate;
    decodeDLSCH.LDPCDecodingAlgorithm = PDSCHExtension.LDPCDecodingAlgorithm;
    decodeDLSCH.MaximumLDPCIterationCount = PDSCHExtension.MaximumLDPCIterationCount;
end

function codedTrBlocks = getDLSCHCodeword(encodeDLSCH,trBlkSizes,Modulation,NumLayers,G,harqEntity)
% Get DL-SCH codeword

    % HARQ processing
    for cwIdx = 1:length(G)
        % If new data for current process and codeword then create a new DL-SCH transport block
        if harqEntity.NewData(cwIdx) 
            trBlk = randi([0 1],trBlkSizes(cwIdx),1);
            setTransportBlock(encodeDLSCH,trBlk,cwIdx-1,harqEntity.HARQProcessID);
        end
    end

    % Encode the DL-SCH transport blocks
    if ~any(Modulation == "4096QAM") 
        codedTrBlocks = encodeDLSCH(Modulation,NumLayers,G,harqEntity.RedundancyVersion,harqEntity.HARQProcessID);
    else
        codedTrBlocks = encodeDLSCH('64QAM',NumLayers,G,harqEntity.RedundancyVersion,harqEntity.HARQProcessID);
        dlschInfo = nrDLSCHInfo(trBlkSizes,encodeDLSCH.TargetCodeRate);            
        codedTrBlocks = encodeReshape(codedTrBlocks,dlschInfo.C,Modulation,NumLayers);
    end

end

function perfEstInfo = perfectEstimatorConfig(pathGains,sampleTimes,pathFilters,noiseEst,perfChEst)
% Perfect channel estimator configuration

    perfEstInfo.PathGains = pathGains;
    perfEstInfo.PathFilters = pathFilters;
    perfEstInfo.SampleTimes = sampleTimes;
    perfEstInfo.NoiseEstimate = noiseEst;
    perfEstInfo.PerfectChannelEstimator = perfChEst;

end

function estChannelGrid = getInitialChannelEstimate(carrier,nTxAnts,propchannel)
% Obtain channel estimate before first transmission. This can be used to
% obtain a precoding matrix for the first slot.

    ofdmInfo = hpre6GOFDMInfo(carrier);
    
    chInfo = info(propchannel);
    maxChDelay = chInfo.MaximumChannelDelay;
    
    % Temporary waveform (only needed for the sizes)
    tmpWaveform = zeros((ofdmInfo.SampleRate/1000/carrier.SlotsPerSubframe)+maxChDelay,nTxAnts,"single");
    
    % Filter through channel    
    [~,pathGains,sampleTimes] = propchannel(tmpWaveform);
    
    % Perfect timing synch    
    pathFilters = getPathFilters(propchannel);
    offset = nrPerfectTimingEstimate(pathGains,pathFilters);
    
    % Perfect channel estimate
    estChannelGrid = hpre6GPerfectChannelEstimate(carrier,pathGains,pathFilters,offset,sampleTimes);
    
end

function plotLayerEVM(NSlots,nslot,pdsch,siz,pdschIndices,pdschSymbols,pdschEq,SNRdB)
% Plot EVM information

    persistent slotEVM;
    persistent rbEVM;
    persistent evmPerSlot;
    
    if (nslot==0)
        slotEVM = comm.EVM;
        rbEVM = comm.EVM;
        evmPerSlot = NaN(NSlots,pdsch.NumLayers);
        figure;
    end
    evmPerSlot(nslot+1,:) = slotEVM(pdschSymbols,pdschEq);
    subplot(2,1,1);
    plot(0:(NSlots-1),evmPerSlot,'o-');
    xlabel("Slot number");
    ylabel("EVM (%)");
    legend("layer " + (1:pdsch.NumLayers),'Location','EastOutside');
    title("EVM per layer per slot. SNR = "+SNRdB+" dB");

    subplot(2,1,2);
    [k,~,p] = ind2sub(siz,pdschIndices);
    rbsubs = floor((k-1) / 12);
    NRB = siz(1) / 12;
    evmPerRB = NaN(NRB,pdsch.NumLayers);
    for nu = 1:pdsch.NumLayers
        for rb = unique(rbsubs).'
            this = (rbsubs==rb & p==nu);
            evmPerRB(rb+1,nu) = rbEVM(pdschSymbols(this),pdschEq(this));
        end
    end
    plot(0:(NRB-1),evmPerRB,'x-');
    xlabel("Resource block");
    ylabel("EVM (%)");
    legend("layer " + (1:pdsch.NumLayers),'Location','EastOutside');
    title("EVM per layer per resource block, slot #"+num2str(nslot)+". SNR = "+SNRdB+" dB");
    
    drawnow;
    
end

function validateParameters(simParameters)

    % Validate the number of layers, relative to the antenna geometry
    numlayers = simParameters.PDSCH.NumLayers;
    ntxants = simParameters.NTxAnts;
    nrxants = simParameters.NRxAnts;
    antennaDescription = sprintf('min(NTxAnts,NRxAnts) = min(%d,%d) = %d',ntxants,nrxants,min(ntxants,nrxants));
    if numlayers > min(ntxants,nrxants)
        error('The number of layers (%d) must satisfy NumLayers <= %s', ...
            numlayers,antennaDescription);
    end
    
    % Display a warning if the maximum possible rank of the channel equals
    % the number of layers
    if (numlayers > 2) && (numlayers == min(ntxants,nrxants))
        warning(['The maximum possible rank of the channel, given by %s, is equal to NumLayers (%d).' ...
            ' This may result in a decoding failure under some channel conditions.' ...
            ' Try decreasing the number of layers or increasing the channel rank' ...
            ' (use more transmit or receive antennas).'],antennaDescription,numlayers); %#ok<SPWRN>
    end

    % Validate modulation schemes for multicodeword transmissions
    modulation = simParameters.PDSCH.Modulation;
    if iscell(modulation) && any(strcmp(modulation,'4096QAM'))
        error("The modulation scheme '4096QAM' is not supported for multicodeword transmissions.")
    end

end

function out = encodeReshape(in,C,Modulation,NumLayers)
    % Undo rate matching for 64QAM
    out = rateReshape("RateRecover",in,C,nr5g.internal.getQm('64QAM'),NumLayers); 

    % Do rate matching for 4096QAM
    out = rateReshape("RateMatch",out,C,nr5g.internal.getQm(Modulation),NumLayers);
end

function out = decodeReshape(in,C,Modulation,NumLayers)
    % Do rate recovery for 4096QAM
    in = in{:};
    in = rateReshape("RateRecover",in,C,nr5g.internal.getQm(Modulation),NumLayers);

    % Undo rate recovery for 64QAM
    out = rateReshape("RateMatch",in,C,nr5g.internal.getQm('64QAM'),NumLayers);
    out = {out};
end

function out = rateReshape(rateMode,in,C,Qm,nlayers)
   
    out = [];
    outlen = length(in);
    idx = 1;

    for r = 1:C 
        if r <= C-mod(outlen/(nlayers*Qm),C)
            E = nlayers*Qm*floor(outlen/(nlayers*Qm*C));
        else
            E = nlayers*Qm*ceil(outlen/(nlayers*Qm*C));
        end

        if rateMode == "RateRecover"
            tmpOut = bitDeinterleaving(in(idx:E+idx-1),E,Qm);
        elseif rateMode == "RateMatch"
            tmpOut = bitInterleaving(in(idx:E+idx-1),E,Qm);
        end
        out = [out;tmpOut];
        idx = idx + E;

    end
 
end

function e = bitInterleaving(e,E,Qm)
    e = reshape(e,E/Qm,Qm);
    e = e.';
    e = e(:);
end

function e = bitDeinterleaving(e,E,Qm)
    e = reshape(e,Qm,E/Qm);
    e = e.';
    e = e(:);
end

function out = getTBS(mod,nlayers,nPRB,NREPerPRB,tcr)
    if ~any(mod == "4096QAM")
        out = nrTBS(mod,nlayers,nPRB,NREPerPRB,tcr);
    else
        out = 2*nrTBS("64QAM",nlayers,nPRB,NREPerPRB,tcr);
    end    
end

function out = getTransportBlockCRC(decodeDLSCH,dlschLLRs,pdsch,harqEntity)

    mod = pdsch.Modulation;
    nlayers = pdsch.NumLayers;
    rv = harqEntity.RedundancyVersion;
    harqID = harqEntity.HARQProcessID;

    if ~any(mod == "4096QAM")
        [~,out] = decodeDLSCH(dlschLLRs,mod,nlayers,rv,harqID);
    else
        dlschInfo = nrDLSCHInfo(decodeDLSCH.TransportBlockLength,decodeDLSCH.TargetCodeRate);
        dlschLLRs = decodeReshape(dlschLLRs,dlschInfo.C,mod,nlayers);
        [~,out] = decodeDLSCH(dlschLLRs,"64QAM",nlayers,rv,harqID);
    end

end