simParameters = [];             % Clear simParameters variable
simParameters.NFrames = 2;      % Number of 10ms frames
simParameters.SNRIn = [-5 0 5]; % SNR range (dB)

perfectChannelEstimator = true;

% Set waveform type and PDSCH numerology (SCS and CP type)
simParameters.NRB = 51;                  % Bandwidth in number of resource blocks (51RBs at 30kHz SCS for 20MHz BW)
simParameters.SubcarrierSpacing = 30;    % 15, 30, 60, 120, 240 (kHz)
simParameters.CyclicPrefix = 'Normal';   % 'Normal' or 'Extended'
simParameters.NCellID = 1;               % Cell identity

% DL-SCH/PDSCH parameters
simParameters.PDSCH.PRBSet = 0:simParameters.NRB-1; % PDSCH PRB allocation
simParameters.PDSCH.SymbolSet = 0:13;           % PDSCH symbol allocation in each slot
simParameters.PDSCH.EnableHARQ = true;          % Enable/disable HARQ, if disabled, single transmission with RV=0, i.e. no retransmissions

simParameters.PDSCH.NLayers = 2;                % Number of PDSCH layers
simParameters.NTxAnts = 8;                      % Number of PDSCH transmission antennas (1,2,4,8,16,32,64,128,256,512,1024) >= NLayers
if simParameters.PDSCH.NLayers > 4              % Multicodeword transmission
    simParameters.NumCW = 2;                        % Number of codewords
    simParameters.PDSCH.TargetCodeRate = [490 490]./1024; % Code rate used to calculate transport block sizes
    simParameters.PDSCH.Modulation = {'16QAM','16QAM'};   % 'QPSK', '16QAM', '64QAM', '256QAM'
    simParameters.NRxAnts = 8;                            % Number of UE receive antennas (even number >= NLayers)
else
    simParameters.NumCW = 1;                        % Number of codewords
    simParameters.PDSCH.TargetCodeRate = 490/1024;  % Code rate used to calculate transport block sizes
    simParameters.PDSCH.Modulation = '16QAM';       % 'QPSK', '16QAM', '64QAM', '256QAM'
    simParameters.NRxAnts = 2;                      % Number of UE receive antennas (1 or even number >= NLayers)
end

% DM-RS and antenna port configuration (TS 38.211 Section 7.4.1.1)
simParameters.PDSCH.PortSet = 0:simParameters.PDSCH.NLayers-1; % DM-RS ports to use for the layers
simParameters.PDSCH.PDSCHMappingType = 'A';     % PDSCH mapping type ('A'(slot-wise),'B'(non slot-wise))
simParameters.PDSCH.DMRSTypeAPosition = 2;      % Mapping type A only. First DM-RS symbol position (2,3)
simParameters.PDSCH.DMRSLength = 1;             % Number of front-loaded DM-RS symbols (1(single symbol),2(double symbol))
simParameters.PDSCH.DMRSAdditionalPosition = 0; % Additional DM-RS symbol positions (max range 0...3)
simParameters.PDSCH.DMRSConfigurationType = 2;  % DM-RS configuration type (1,2)
simParameters.PDSCH.NumCDMGroupsWithoutData = 1;% Number of CDM groups without data
simParameters.PDSCH.NIDNSCID = 1;               % Scrambling identity (0...65535)
simParameters.PDSCH.NSCID = 0;                  % Scrambling initialization (0,1)
% Reserved PRB patterns (for CORESETs, forward compatibility etc)
simParameters.PDSCH.Reserved.Symbols = [];      % Reserved PDSCH symbols
simParameters.PDSCH.Reserved.PRB = [];          % Reserved PDSCH PRBs
simParameters.PDSCH.Reserved.Period = [];       % Periodicity of reserved resources
% PDSCH resource block mapping (TS 38.211 Section 7.3.1.6)
simParameters.PDSCH.VRBToPRBInterleaving = 0;   % Disable interleaved resource mapping
simParameters.PDSCH.VRBBundleSize = 4;
% PDSCH PRB bundling (TS 38.214 Section 5.1.2.3)
simParameters.PDSCH.PRGBundleSize = [];         % 2, 4, or [] to signify "wideband"
% PT-RS configuration (TS 38.211 Section 7.4.1.2)
simParameters.PDSCH.EnablePTRS = 0;             % Enable or disable PT-RS (1 or 0)
simParameters.PDSCH.PTRSTimeDensity = 1;        % PT-RS time density (L_PT-RS) (1, 2, 4)
simParameters.PDSCH.PTRSFrequencyDensity = 2;   % PT-RS frequency density (K_PT-RS) (2 or 4)
simParameters.PDSCH.PTRSREOffset = '00';        % PT-RS resource element offset ('00', '01', '10', '11')
simParameters.PDSCH.PTRSPortSet = [];           % PT-RS antenna port, subset of DM-RS port set. Empty corresponds to lower DM-RS port number

% Define the propagation channel type
simParameters.ChannelType = 'CDL'; % 'CDL' or 'TDL'

% SS burst configuration
simParameters.SSBurst.BlockPattern = 'Case B';    % 30kHz subcarrier spacing
simParameters.SSBurst.SSBTransmitted = [0 1 0 1]; % Bitmap indicating blocks transmitted in the burst
simParameters.SSBurst.SSBPeriodicity = 20;        % SS burst set periodicity in ms (5, 10, 20, 40, 80, 160)

validateNLayers(simParameters);

carrier = nrCarrierConfig;
carrier.SubcarrierSpacing = simParameters.SubcarrierSpacing;
carrier.CyclicPrefix = simParameters.CyclicPrefix;
carrier.NSizeGrid = simParameters.NRB;
carrier.NCellID = simParameters.NCellID;
pdsch = simParameters.PDSCH;

% Specify additional required fields for PDSCH
pdsch.RNTI = 1;

Xoh_PDSCH = 6*simParameters.PDSCH.EnablePTRS;     % Set Xoh-PDSCH overhead to 0 when PT-RS is disabled or to 6 otherwise

snrIn = simParameters.SNRIn;
nTxAnts = simParameters.NTxAnts;
nRxAnts = simParameters.NRxAnts;
channelType = simParameters.ChannelType;

if strcmpi(channelType,'CDL')
    channel = nrCDLChannel; % CDL channel object

    % Use CDL-C model (Urban macrocell model)
    channel.DelayProfile = 'CDL-C';
    channel.DelaySpread = 300e-9;

    % Turn the overall number of antennas into a specific antenna panel
    % array geometry. The number of antennas configured is updated when
    % nTxAnts is not one of (1,2,4,8,16,32,64,128,256,512,1024) or nRxAnts
    % is not 1 or even.
    [channel.TransmitAntennaArray.Size, channel.ReceiveAntennaArray.Size] = ...
        hArrayGeometry(nTxAnts,nRxAnts);
    nTxAnts = prod(channel.TransmitAntennaArray.Size);
    nRxAnts = prod(channel.ReceiveAntennaArray.Size);
    simParameters.NTxAnts = nTxAnts;
    simParameters.NRxAnts = nRxAnts;

elseif strcmpi(channelType,'TDL')
    channel = nrTDLChannel; % TDL channel object
    % Set the channel geometry
    channel.DelayProfile = 'TDL-C';
    channel.DelaySpread = 300e-9;
    channel.NumTransmitAntennas = nTxAnts;
    channel.NumReceiveAntennas = nRxAnts;
else
    error('ChannelType parameter field must be either CDL or TDL.');
end

waveformInfo = nrOFDMInfo(carrier);
channel.SampleRate = waveformInfo.SampleRate;

chInfo = info(channel);
maxChDelay = ceil(max(chInfo.PathDelays*channel.SampleRate)) + chInfo.ChannelFilterDelay;

% Create SS burst information structure
ssburst = simParameters.SSBurst;
ssburst.NCellID = carrier.NCellID;
ssburst.SampleRate = waveformInfo.SampleRate;

ssbInfo = hSSBurstInfo(ssburst);

% Map the occupied subcarriers and transmitted symbols of the SS burst
% (defined in the SS burst numerology) to PDSCH PRBs and symbols in the
% data numerology
[mappedPRB,mappedSymbols] = mapNumerology(ssbInfo.OccupiedSubcarriers,ssbInfo.OccupiedSymbols,ssbInfo.NRB,carrier.NSizeGrid,ssbInfo.SubcarrierSpacing,carrier.SubcarrierSpacing);

% Configure the PDSCH to reserve these resources so that the PDSCH
% transmission does not overlap the SS burst
reservation.Symbols = mappedSymbols;
reservation.PRB = mappedPRB;
reservation.Period = ssburst.SSBPeriodicity * (carrier.SubcarrierSpacing/15); % Period in slots
pdsch.Reserved(end+1) = reservation;

% Array to store the maximum throughput for all SNR points
maxThroughput = zeros(length(snrIn),1);
% Array to store the simulation throughput for all SNR points
simThroughput = zeros(length(snrIn),1);

% Set up Redundancy Version (RV) sequence, number of HARQ processes and
% the sequence in which the HARQ processes are used
if pdsch.EnableHARQ
    % In the final report of RAN WG1 meeting #91 (R1-1719301), it was
    % observed in R1-1717405 that if performance is the priority, [0 2 3 1]
    % should be used. If self-decodability is the priority, it should be
    % taken into account that the upper limit of the code rate at which
    % each RV is self-decodable is in the following order: 0>3>>2>1
    rvSeq = [0 2 3 1];
else
    % HARQ disabled - single transmission with RV=0, no retransmissions
    rvSeq = 0;
end

% Create DLSCH encoder system object
encodeDLSCH = nrDLSCH;
encodeDLSCH.MultipleHARQProcesses = true;
encodeDLSCH.TargetCodeRate = pdsch.TargetCodeRate;

% Create DLSCH decoder system object
% Use layered belief propagation for LDPC decoding, with half the number of
% iterations as compared to the default for belief propagation decoding
decodeDLSCH = nrDLSCHDecoder;
decodeDLSCH.MultipleHARQProcesses = true;
decodeDLSCH.TargetCodeRate = pdsch.TargetCodeRate;
decodeDLSCH.LDPCDecodingAlgorithm = 'Layered belief propagation';
decodeDLSCH.MaximumLDPCIterationCount = 6;

% The temporary variables 'carrier_init', 'pdsch_init', 'ssburst_init' and
% 'decodeDLSCH_init' are used to create the temporary variables 'carrier',
% 'pdsch', 'ssburst' and 'decodeDLSCH' within the SNR loop to create
% independent instances in case of parallel simulation
carrier_init = carrier;
pdsch_init = pdsch;
ssburst_init = ssburst;
decodeDLSCH_init = clone(decodeDLSCH);

% Get the values of 'NFrames' and 'NumCW' to avoid overhead in case of
% parallel simulation
NFrames = simParameters.NFrames;
NumCW = simParameters.NumCW;

for snrIdx = 1:numel(snrIn) % comment out for parallel computing
% parfor snrIdx = 1:numel(snrIn) % uncomment for parallel computing
% To reduce the total simulation time, you can execute this loop in
% parallel by using the Parallel Computing Toolbox. Comment out the 'for'
% statement and uncomment the 'parfor' statement. If the Parallel Computing
% Toolbox is not installed, 'parfor' defaults to normal 'for' statement

    % Set the random number generator settings to default values
    rng('default');

    % Initialize variables for this SNR point, required for initialization
    % of variables when using the Parallel Computing Toolbox
    carrier = carrier_init;
    pdsch = pdsch_init;
    ssburst = ssburst_init;
    decodeDLSCH = clone(decodeDLSCH_init);
    pathFilters = [];
    ssbWaveform = [];

    SNRdB = snrIn(snrIdx);
    fprintf('\nSimulating transmission scheme 1 (%dx%d) and SCS=%dkHz with %s channel at %gdB SNR for %d 10ms frame(s)\n',...
        nTxAnts,nRxAnts,carrier.SubcarrierSpacing, ...
        channelType,SNRdB,NFrames);

    % Initialize variables used in the simulation and analysis
    bitTput = [];           % Number of successfully received bits per transmission
    txedTrBlkSizes = [];    % Number of transmitted info bits per transmission

    % Specify the order in which we cycle through the HARQ processes
    NHARQProcesses = 16;
    harqSequence = 1:NHARQProcesses;

    % Initialize the state of all HARQ processes
    harqProcesses = hNewHARQProcesses(NHARQProcesses,rvSeq,NumCW);
    harqProcCntr = 0; % HARQ process counter

    % Reset the channel so that each SNR point will experience the same
    % channel realization
    reset(channel);

    % Total number of slots in the simulation period
    NSlots = NFrames * carrier.SlotsPerFrame;

    % Create 'gnb' structure
    gnb = struct();
    gnb.NRB = carrier.NSizeGrid;
    gnb.CyclicPrefix = carrier.CyclicPrefix;
    gnb.SubcarrierSpacing = carrier.SubcarrierSpacing;

    % Index to the start of the current set of SS burst samples to be
    % transmitted
    ssbSampleIndex = 1;

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

    % Timing offset, updated in every slot for perfect synchronization and
    % when the correlation is strong for practical synchronization
    offset = 0;

    % Loop over the entire waveform length
    for nslot = 0:NSlots-1

        % Update the carrier and PDSCH slot numbers to account for a new
        % PDSCH transmission
        carrier.NSlot = nslot;
        pdsch.NSlot = nslot;

        % Generate a new SS burst when necessary
        if (ssbSampleIndex==1)
            nSubframe = carrier.NSlot / carrier.SlotsPerSubframe;
            ssburst.NFrame = floor(nSubframe / 10);
            ssburst.NHalfFrame = mod(nSubframe / 5,2);
            [ssbWaveform,~,ssbInfo] = hSSBurst(ssburst);
        end

        % Get HARQ process index for the current PDSCH from HARQ index table
        harqProcIdx = harqSequence(mod(harqProcCntr,length(harqSequence))+1);

        % Update current HARQ process information (this updates the RV
        % depending on CRC pass or fail in the previous transmission for
        % this HARQ process)
        harqProcesses(harqProcIdx) = hUpdateHARQProcess(harqProcesses(harqProcIdx),NumCW);

        % Calculate the transport block sizes for the codewords in the slot
        [pdschIndices,dmrsIndices,dmrsSymbols,ptrsIndices,ptrsSymbols,pdschIndicesInfo] = hPDSCHResources(gnb,pdsch);
        trBlkSizes = nrTBS(pdsch.Modulation,pdsch.NLayers,numel(pdsch.PRBSet),pdschIndicesInfo.NREPerPRB,pdsch.TargetCodeRate,Xoh_PDSCH);

        % HARQ: check CRC from previous transmission per codeword, i.e. is
        % a retransmission required?
        for cwIdx = 1:NumCW
            NDI = false;
            if harqProcesses(harqProcIdx).blkerr(cwIdx) % Errored
                if (harqProcesses(harqProcIdx).RVIdx(cwIdx)==1) % end of rvSeq
                    resetSoftBuffer(decodeDLSCH,cwIdx-1,harqProcIdx-1);
                    NDI = true;
                end
            else    % No error
                NDI = true;
            end
            if NDI
                trBlk = randi([0 1],trBlkSizes(cwIdx),1);
                setTransportBlock(encodeDLSCH,trBlk,cwIdx-1,harqProcIdx-1);
            end
        end

        % Encode the DL-SCH transport blocks
        codedTrBlock = encodeDLSCH(pdsch.Modulation,pdsch.NLayers,...
            pdschIndicesInfo.G,harqProcesses(harqProcIdx).RV,harqProcIdx-1);

        % Get wtx (precoding matrix) calculated in previous slot
        wtx = newWtx;

        % PDSCH modulation and precoding
        pdschSymbols = nrPDSCH(codedTrBlock,pdsch.Modulation,pdsch.NLayers,carrier.NCellID,pdsch.RNTI);
        pdschGrid = nrResourceGrid(carrier,nTxAnts);
        [pdschSymbols,pdschAntIndices] = hPRGPrecode(size(pdschGrid),carrier.NStartGrid,pdschSymbols,pdschIndices,wtx);

        % PDSCH mapping in grid associated with PDSCH transmission period
        pdschGrid(pdschAntIndices) = pdschSymbols;

        % PDSCH DM-RS precoding and mapping
        [dmrsAntSymbols,dmrsAntIndices] = hPRGPrecode(size(pdschGrid),carrier.NStartGrid,dmrsSymbols,dmrsIndices,wtx);
        pdschGrid(dmrsAntIndices) = dmrsAntSymbols;

        % PDSCH PT-RS precoding and mapping
        [ptrsAntSymbols,ptrsAntIndices] = hPRGPrecode(size(pdschGrid),carrier.NStartGrid,ptrsSymbols,ptrsIndices,wtx);
        pdschGrid(ptrsAntIndices) = ptrsAntSymbols;

        % OFDM modulation of associated resource elements
        txWaveform = nrOFDMModulate(carrier, pdschGrid);

        % Add the appropriate portion of SS burst waveform to the
        % transmitted waveform
        Nt = size(txWaveform,1);
        txWaveform = txWaveform + ssbWaveform(ssbSampleIndex + (0:Nt-1),:);
        ssbSampleIndex = mod(ssbSampleIndex + Nt,size(ssbWaveform,1));

        % Pass data through channel model. Append zeros at the end of the
        % transmitted waveform to flush channel content. These zeros take
        % into account any delay introduced in the channel. This is a mix
        % of multipath delay and implementation delay. This value may
        % change depending on the sampling rate, delay profile and delay
        % spread
        txWaveform = [txWaveform; zeros(maxChDelay, size(txWaveform,2))];
        [rxWaveform,pathGains,sampleTimes] = channel(txWaveform);

        % 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/20); % Calculate linear noise gain
        N0 = 1/(sqrt(2.0*nRxAnts*double(waveformInfo.Nfft))*SNR);
        noise = N0*complex(randn(size(rxWaveform)),randn(size(rxWaveform)));
        rxWaveform = rxWaveform + noise;

        if (perfectChannelEstimator)
            % Perfect synchronization. Use information provided by the channel
            % to find the strongest multipath component
            pathFilters = getPathFilters(channel); % get path filters for perfect channel estimation
            [offset,mag] = nrPerfectTimingEstimate(pathGains,pathFilters);
        else
            % Practical synchronization. Correlate the received waveform
            % with the PDSCH DM-RS to give timing offset estimate 't' and
            % correlation magnitude 'mag'. The function
            % hSkipWeakTimingOffset is used to update the receiver timing
            % offset. If the correlation peak in 'mag' is weak, the current
            % timing estimate 't' is ignored and the previous estimate
            % 'offset' is used
            [t,mag] = nrTimingEstimate(carrier,rxWaveform,dmrsIndices,dmrsSymbols); %#ok<UNRCH>
            offset = hSkipWeakTimingOffset(offset,t,mag);
        end
        rxWaveform = rxWaveform(1+offset:end, :);

        % Perform OFDM demodulation on the received data to recreate the
        % resource grid, including padding in the event that practical
        % synchronization results in an incomplete slot being demodulated
        rxGrid = nrOFDMDemodulate(carrier, rxWaveform);
        [K,L,R] = size(rxGrid);
        if (L < carrier.SymbolsPerSlot)
            rxGrid = cat(2,rxGrid,zeros(K,carrier.SymbolsPerSlot-L,R));
        end

        if (perfectChannelEstimator)
            % Perfect channel estimation, using the value of the path gains
            % provided by the channel. This channel estimate does not
            % include the effect of transmitter precoding
            estChannelGrid = nrPerfectChannelEstimate(carrier,pathGains,pathFilters,offset,sampleTimes);

            % Get perfect noise estimate (from the noise realization)
            noiseGrid = nrOFDMDemodulate(carrier,noise(1+offset:end ,:));
            noiseEst = var(noiseGrid(:));

            % Get precoding matrix for next slot
            newWtx = getPrecodingMatrix(carrier,pdsch,estChannelGrid);

            % Get PDSCH resource elements from the received grid and
            % channel estimate
            [pdschRx,pdschHest,~,pdschHestIndices] = nrExtractResources(pdschIndices,rxGrid,estChannelGrid);

            % Apply precoding to channel estimate
            pdschHest = hPRGPrecode(size(estChannelGrid),carrier.NStartGrid,pdschHest,pdschHestIndices,permute(wtx,[2 1 3]));
        else
            % Practical channel estimation between the received grid and
            % each transmission layer, using the PDSCH DM-RS for each
            % layer. This channel estimate includes the effect of
            % transmitter precoding
            [estChannelGrid,noiseEst] = nrChannelEstimate(carrier,rxGrid,dmrsIndices,dmrsSymbols,'CDMLengths',pdschIndicesInfo.CDMLengths); %#ok<UNRCH>

            % Get PDSCH resource elements from the received grid and
            % channel estimate
            [pdschRx,pdschHest] = nrExtractResources(pdschIndices,rxGrid,estChannelGrid);

            % Remove precoding from estChannelGrid prior to precoding
            % matrix calculation
            estChannelGridPorts = precodeChannelEstimate(carrier,estChannelGrid,conj(wtx));

            % Get precoding matrix for next slot
            newWtx = getPrecodingMatrix(carrier,pdsch,estChannelGridPorts);
        end

        % Equalization
        [pdschEq,csi] = nrEqualizeMMSE(pdschRx,pdschHest,noiseEst);

        % Common phase error (CPE) compensation
        if ~isempty(ptrsIndices)
            % Initialize temporary grid to store equalized symbols
            tempGrid = nrResourceGrid(carrier,pdsch.NLayers);

            % Extract PT-RS symbols from received grid and estimated
            % channel grid
            [ptrsRx,ptrsHest,~,~,ptrsHestIndices,ptrsLayerIndices] = nrExtractResources(ptrsIndices,rxGrid,estChannelGrid,tempGrid);

            if (perfectChannelEstimator)
                % Apply precoding to channel estimate
                ptrsHest = hPRGPrecode(size(estChannelGrid),carrier.NStartGrid,ptrsHest,ptrsHestIndices,permute(wtx,[2 1 3]));
            end

            % Equalize PT-RS symbols and map them to tempGrid
            ptrsEq = nrEqualizeMMSE(ptrsRx,ptrsHest,noiseEst);
            tempGrid(ptrsLayerIndices) = ptrsEq;

            % Estimate the residual channel at the PT-RS locations in
            % tempGrid
            cpe = nrChannelEstimate(tempGrid,ptrsIndices,ptrsSymbols);

            % Sum estimates across subcarriers, receive antennas, and
            % layers. Then, get the CPE by taking the angle of the
            % resultant sum
            cpe = angle(sum(cpe,[1 3 4]));

            % Map the equalized PDSCH symbols to tempGrid
            tempGrid(pdschIndices) = pdschEq;

            % Correct CPE in each OFDM symbol within the range of reference
            % PT-RS OFDM symbols
            symLoc = pdschIndicesInfo.PTRSSymbolSet(1)+1:pdschIndicesInfo.PTRSSymbolSet(end)+1;
            tempGrid(:,symLoc,:) = tempGrid(:,symLoc,:).*exp(-1i*cpe(symLoc));

            % Extract PDSCH symbols
            pdschEq = tempGrid(pdschIndices);
        end

        % Decode PDSCH physical channel
        [dlschLLRs,rxSymbols] = nrPDSCHDecode(pdschEq,pdsch.Modulation,carrier.NCellID,pdsch.RNTI,noiseEst);

        % Scale LLRs by CSI
        csi = nrLayerDemap(csi); % CSI layer demapping
        for cwIdx = 1:NumCW
            Qm = length(dlschLLRs{cwIdx})/length(rxSymbols{cwIdx}); % bits per symbol
            csi{cwIdx} = repmat(csi{cwIdx}.',Qm,1);   % expand by each bit per symbol
            dlschLLRs{cwIdx} = dlschLLRs{cwIdx} .* csi{cwIdx}(:);   % scale
        end

        % Decode the DL-SCH transport channel
        decodeDLSCH.TransportBlockLength = trBlkSizes;
        [decbits,harqProcesses(harqProcIdx).blkerr] = decodeDLSCH(dlschLLRs,pdsch.Modulation,pdsch.NLayers,harqProcesses(harqProcIdx).RV,harqProcIdx-1);

        % Store values to calculate throughput (only for active PDSCH instances)
        if(any(trBlkSizes ~= 0))
            bitTput = [bitTput trBlkSizes.*(1-harqProcesses(harqProcIdx).blkerr)];
            txedTrBlkSizes = [txedTrBlkSizes trBlkSizes];
        end

        % Update HARQ process counter
        harqProcCntr = harqProcCntr + 1;

        % Display transport block error information per codeword managed by current HARQ process
        fprintf('\n(%3.2f%%) HARQ Proc %d: ',100*(nslot+1)/NSlots,harqProcIdx);
        estrings = {'passed','failed'};
        rvi = harqProcesses(harqProcIdx).RVIdx;
        for cw=1:length(rvi)
            cwrvi = rvi(cw);
            % Create a report on the RV state given position in RV sequence and decoding error
            if cwrvi == 1
                ts = sprintf('Initial transmission (RV=%d)',rvSeq(cwrvi));
            else
                ts = sprintf('Retransmission #%d (RV=%d)',cwrvi-1,rvSeq(cwrvi));
            end
            fprintf('CW%d:%s %s. ',cw-1,ts,estrings{1+harqProcesses(harqProcIdx).blkerr(cw)});
        end

     end

    % Calculate maximum and simulated throughput
    maxThroughput(snrIdx) = sum(txedTrBlkSizes); % Max possible throughput
    simThroughput(snrIdx) = sum(bitTput,2);      % Simulated throughput

    % Display the results dynamically in the command window
    fprintf([['\n\nThroughput(Mbps) for ', num2str(NFrames) ' frame(s) '],...
        '= %.4f\n'], 1e-6*simThroughput(snrIdx)/(NFrames*10e-3));
    fprintf(['Throughput(%%) for ', num2str(NFrames) ' frame(s) = %.4f\n'],...
        simThroughput(snrIdx)*100/maxThroughput(snrIdx));

end

figure;
plot(snrIn,simThroughput*100./maxThroughput,'o-.')
xlabel('SNR (dB)'); ylabel('Throughput (%)'); grid on;
title(sprintf('(%dx%d) / NRB=%d / SCS=%dkHz',...
              nTxAnts,nRxAnts,carrier_init.NSizeGrid,carrier_init.SubcarrierSpacing));

% Bundle key parameters and results into a combined structure for recording
simResults.simParameters = simParameters;
simResults.simThroughput = simThroughput;

function validateNLayers(simParameters)
% Validate the number of layers
    if length(simParameters.PDSCH.PortSet)~= simParameters.PDSCH.NLayers
        error('The number of elements of PortSet (%d) must be the same as the number of layers (%d)',...
            length(simParameters.PDSCH.PortSet), simParameters.PDSCH.NLayers);
    end

    if simParameters.PDSCH.NLayers > min(simParameters.NTxAnts,simParameters.NRxAnts)
        error('The number of layers (%d) must satisfy NLayers <= min(NTxAnts,NRxAnts) = min(%d,%d) = (%d)',...
            simParameters.PDSCH.NLayers,simParameters.NTxAnts,simParameters.NRxAnts,min(simParameters.NTxAnts,simParameters.NRxAnts));
    end
end

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

    ofdmInfo = nrOFDMInfo(carrier);

    chInfo = info(channel);
    maxChDelay = ceil(max(chInfo.PathDelays*channel.SampleRate)) + chInfo.ChannelFilterDelay;

    % Temporary waveform (only needed for the sizes)
    tmpWaveform = zeros((ofdmInfo.SampleRate/1000/carrier.SlotsPerSubframe)+maxChDelay,nTxAnts);

    % Filter through channel
    [~,pathGains,sampleTimes] = channel(tmpWaveform);

    % Perfect timing synch
    pathFilters = getPathFilters(channel);
    offset = nrPerfectTimingEstimate(pathGains,pathFilters);

    % Perfect channel estimate
    estChannelGrid = nrPerfectChannelEstimate(carrier,pathGains,pathFilters,offset,sampleTimes);

end

function wtx = getPrecodingMatrix(carrier,pdsch,hestGrid)
% Calculate precoding matrices for all PRGs in the carrier that overlap
% with the PDSCH allocation

    % Maximum CRB addressed by carrier grid
    maxCRB = carrier.NStartGrid + carrier.NSizeGrid - 1;

    % PRG size
    if (isfield(pdsch,'PRGBundleSize') && ~isempty(pdsch.PRGBundleSize))
        Pd_BWP = pdsch.PRGBundleSize;
    else
        Pd_BWP = maxCRB + 1;
    end

    % PRG numbers (1-based) for each RB in the carrier grid
    NPRG = ceil((maxCRB + 1) / Pd_BWP);
    prgset = repmat((1:NPRG),Pd_BWP,1);
    prgset = prgset(carrier.NStartGrid + (1:carrier.NSizeGrid).');

    [~,~,R,P] = size(hestGrid);
    wtx = zeros([pdsch.NLayers P NPRG]);
    for i = 1:NPRG

        % Subcarrier indices within current PRG and within the PDSCH
        % allocation
        thisPRG = find(prgset==i) - 1;
        thisPRG = intersect(thisPRG,pdsch.PRBSet(:) + carrier.NStartGrid,'rows');
        prgSc = (1:12)' + 12*thisPRG';
        prgSc = prgSc(:);

        if (~isempty(prgSc))

            % Average channel estimate in PRG
            estAllocGrid = hestGrid(prgSc,:,:,:);
            Hest = permute(mean(reshape(estAllocGrid,[],R,P)),[2 3 1]);

            % SVD decomposition
            [~,~,V] = svd(Hest);
            wtx(:,:,i) = V(:,1:pdsch.NLayers).';

        end

    end

    wtx = wtx / sqrt(pdsch.NLayers); % Normalize by NLayers

end

function estChannelGrid = precodeChannelEstimate(carrier,estChannelGrid,W)
% Apply precoding matrix W to the last dimension of the channel estimate

    [K,L,R,P] = size(estChannelGrid);
    estChannelGrid = reshape(estChannelGrid,[K*L R P]);
    estChannelGrid = hPRGPrecode([K L R P],carrier.NStartGrid,estChannelGrid,reshape(1:numel(estChannelGrid),[K*L R P]),W);
    estChannelGrid = reshape(estChannelGrid,K,L,R,[]);

end

function [mappedPRB,mappedSymbols] = mapNumerology(subcarriers,symbols,nrbs,nrbt,fs,ft)
% Map the SSBurst numerology to PDSCH numerology. The outputs are:
%   - mappedPRB: 0-based PRB indices for carrier resource grid (arranged in a column)
%   - mappedSymbols: 0-based OFDM symbol indices in a slot for carrier resource grid (arranged in a row)
%     carrier resource grid is sized using gnb.NRB, gnb.CyclicPrefix, spanning 1 slot
% The input parameters are:
%   - subcarriers: 1-based row subscripts for SSB resource grid (arranged in a column)
%   - symbols: 1-based column subscripts for SSB resource grid (arranged in an N-by-4 matrix, 4 symbols for each transmitted burst in a row, N transmitted bursts)
%     SSB resource grid is sized using ssbInfo.NRB, normal CP, spanning 5 subframes
%   - nrbs: source (SSB) NRB
%   - nrbt: target (carrier) NRB
%   - fs: source (SSB) SCS
%   - ft: target (carrier) SCS

    mappedPRB = unique(fix((subcarriers-(nrbs*6) - 1)*fs/(ft*12) + nrbt/2),'stable');

    symbols = symbols.';
    symbols = symbols(:).' - 1;

    if (ft < fs)
        % If ft/fs < 1, reduction
        mappedSymbols = unique(fix(symbols*ft/fs),'stable');
    else
        % Else, repetition by ft/fs
        mappedSymbols = reshape((0:(ft/fs-1))' + symbols(:)'*ft/fs,1,[]);
    end

end