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

displaySimulationInformation = true;

perfectChannelEstimator = true;

% Bandwidth, numerology (SCS and CP type) and other general parameters
simParameters.NRB = 52;                % Bandwidth in number of resource blocks (52RBs at 15kHz SCS for 10MHz BW)
simParameters.SubcarrierSpacing = 15;  % 15, 30, 60, 120, 240 (kHz)
simParameters.CyclicPrefix = 'Normal'; % 'Normal' or 'Extended'
simParameters.NCellID = 0;             % Cell identity
simParameters.NTxAnts = 1;             % Number of transmit antennas
simParameters.NRxAnts = 2;             % Number of receive antennas

% UL-SCH/PUSCH parameters
simParameters.PUSCH.TargetCodeRate = 193 / 1024;      % Code rate used to calculate transport block sizes
simParameters.PUSCH.PRBSet = (0:simParameters.NRB-1); % PUSCH PRB allocation
simParameters.PUSCH.SymbolSet = 0:13;            % PUSCH symbol allocation in each slot
simParameters.PUSCH.NohPRB = 0;                  % Additional RE overhead per PRB
simParameters.PUSCH.EnableHARQ = true;           % Enable/disable HARQ, if disabled, single transmission with RV=0, i.e. no retransmissions
simParameters.PUSCH.Modulation = 'QPSK';         % 'pi/2-BPSK', 'QPSK', '16QAM', '64QAM', '256QAM'
simParameters.PUSCH.NLayers = 1;                 % Number of PUSCH layers
simParameters.PUSCH.RNTI = 1;                    % Radio Network Temporary Identifier
simParameters.PUSCH.TransformPrecoding = false;  % Enable/disable transform precoding
simParameters.PUSCH.TxScheme = 'nonCodebook';    % Transmission scheme ('nonCodebook','codebook')
simParameters.PUSCH.NAntennaPorts = 1;           % Number of antenna ports for codebook based precoding
simParameters.PUSCH.TPMI = 0;                    % Precoding matrix indicator for codebook based precoding
% PUSCH DM-RS configuration
simParameters.PUSCH.PUSCHMappingType = 'A';      % PUSCH mapping type ('A'(slot-wise),'B'(non slot-wise))
simParameters.PUSCH.DMRSTypeAPosition = 2;       % Mapping type A only. First DM-RS symbol position (2,3)
simParameters.PUSCH.DMRSLength = 1;              % Number of front-loaded DM-RS symbols (1(single symbol),2(double symbol))
simParameters.PUSCH.DMRSAdditionalPosition = 1;  % Additional DM-RS symbol positions (max range 0...3)
simParameters.PUSCH.DMRSConfigurationType = 1;   % DM-RS configuration type (1,2)
simParameters.PUSCH.NumCDMGroupsWithoutData = 2; % CDM groups without data
simParameters.PUSCH.NIDNSCID = 0;                % Scrambling identity (0...65535)
simParameters.PUSCH.NSCID = 0;                   % Scrambling initialization (0,1)
simParameters.PUSCH.NRSID = 0;                   % Scrambling ID for low-PAPR sequences (0...1007)
simParameters.PUSCH.GroupHopping = 'Disable';    % Hopping type ('Enable','Disable')

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

ue = simParameters;
pusch = simParameters.PUSCH;

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

if strcmpi(channelType,'CDL')

    channel = nrCDLChannel;
    channel.DelayProfile = 'CDL-A';
    [txsize,rxsize] = hArrayGeometry(nTxAnts,nRxAnts,'uplink');
    channel.TransmitAntennaArray.Size = txsize;
    channel.ReceiveAntennaArray.Size = rxsize;

else

    channel = nrTDLChannel;
    channel.DelayProfile = 'TDL-A';
    channel.NumTransmitAntennas = nTxAnts;
    channel.NumReceiveAntennas = nRxAnts;

end

channel.DelaySpread = 30e-9; % in seconds
channel.MaximumDopplerShift = 10; % in Hz

waveformInfo = hOFDMInfo(ue);
channel.SampleRate = waveformInfo.SamplingRate;

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

% 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 pusch.EnableHARQ
    % From PUSCH demodulation requirements in RAN WG4 meeting #88bis
    % (R4-1814062)
    rvSeq = [0 2 3 1];
else
    % HARQ disabled - single transmission with RV=0, no retransmissions
    rvSeq = 0;
end

% Create UL-SCH encoder System object
encodeULSCH = nrULSCH;
encodeULSCH.MultipleHARQProcesses = true;
encodeULSCH.TargetCodeRate = pusch.TargetCodeRate;

% Create UL-SCH 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
decodeULSCH = nrULSCHDecoder;
decodeULSCH.MultipleHARQProcesses = true;
decodeULSCH.TargetCodeRate = pusch.TargetCodeRate;
decodeULSCH.LDPCDecodingAlgorithm = 'Layered belief propagation';
decodeULSCH.MaximumLDPCIterationCount = 6;

% The temporary variables 'ue_init', 'pusch_init' and 'decodeULSCH_init'
% are used to create the temporary variables 'ue', 'pusch' and
% 'decodeULSCH' within the SNR loop to create independent instances in case
% of parallel simulation
ue_init = ue;
pusch_init = pusch;
decodeULSCH_init = clone(decodeULSCH);

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.
% Because parfor-loop iterations are executed in parallel in a
% nondeterministic order, the simulation information displayed for each SNR
% point can be intertwined. To switch off simulation information display,
% set the 'displaySimulationInformation' variable above to false

    % Reset the random number generator and channel, so that each SNR point
    % will experience the same noise and channel realizations
    rng('default');
    reset(channel);

    % Initialize variables for this SNR point, required for initialization
    % of variables when using the Parallel Computing Toolbox
    ue = ue_init;
    pusch = pusch_init;
    pathFilters = [];

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

    % Initialize the state of all HARQ processes and reset the UL-SCH
    % decoder
    harqProcesses = hNewHARQProcesses(NHARQProcesses,rvSeq,1);
    harqProcCntr = 0; % HARQ process counter
    decodeULSCH = clone(decodeULSCH_init);

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

    % Total number of OFDM symbols in the simulation period
    waveformInfo = hOFDMInfo(ue);
    NSymbols = ue.NFrames * 10 * waveformInfo.SymbolsPerSubframe;

    % OFDM symbol number associated with start of each PUSCH transmission
    ue.NSymbol = 0;

    % Running counter of the number of PUSCH transmission instances
    % The simulation will use this counter as the slot number for each
    % PUSCH
    pusch.NSlot = 0;

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

    while ue.NSymbol < NSymbols

        % Calculate the transport block size for this slot
        [puschIndices,dmrsIndices,dmrsSymbols,puschIndicesInfo] = hPUSCHResources(ue,pusch);
        TBS = hPUSCHTBS(pusch,puschIndicesInfo.NREPerPRB - pusch.NohPRB);

        % Get HARQ process index for the current PUSCH from the 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),1);

        % HARQ: check CRC from previous transmission, i.e. is a
        % retransmission required?
        NDI = false;
        if harqProcesses(harqProcIdx).blkerr % errored
            if (harqProcesses(harqProcIdx).RVIdx==1) % end of rvSeq
                resetSoftBuffer(decodeULSCH,harqProcIdx-1);
                NDI = true;
            end
        else % no error
            NDI = true;
        end
        if NDI
            trBlk = randi([0 1],TBS,1);
            setTransportBlock(encodeULSCH,trBlk,harqProcIdx-1);
        end

        % UL-SCH encoding
        codedTrBlock = encodeULSCH(pusch.Modulation,pusch.NLayers,puschIndicesInfo.G,harqProcesses(harqProcIdx).RV,harqProcIdx-1);

        % PUSCH modulation, including codebook based MIMO precoding if
        % TxScheme = 'codebook'
        MRB = numel(pusch.PRBSet);
        puschSymbols = nrPUSCH(codedTrBlock,pusch.Modulation,pusch.NLayers,ue.NCellID,pusch.RNTI,pusch.TransformPrecoding,MRB,pusch.TxScheme,pusch.NAntennaPorts,pusch.TPMI);

        % Create resource grid associated with PUSCH transmission period
        puschGrid = zeros(waveformInfo.NSubcarriers,waveformInfo.SymbolsPerSlot,nTxAnts);

        % Implementation-specific PUSCH MIMO precoding and mapping. This
        % MIMO precoding step is in addition to any codebook based
        % MIMO precoding done during PUSCH modulation above
        if (strcmpi(pusch.TxScheme,'codebook'))
            % codebook based MIMO precoding, F precodes between PUSCH
            % transmit antenna ports and transmit antennas
            F = eye(pusch.NAntennaPorts,nTxAnts);
        else
            % non-codebook based MIMO precoding, F precodes between PUSCH
            % layers and transmit antennas
            F = eye(pusch.NLayers,nTxAnts);
        end
        [~,puschAntIndices] = nrExtractResources(puschIndices,puschGrid);
        puschGrid(puschAntIndices) = puschSymbols * F;

        % Implementation-specific PUSCH DM-RS MIMO precoding and mapping.
        % The DM-RS creation in hPUSCHResources above includes codebook
        % based MIMO precoding if applicable
        for p = 1:size(dmrsSymbols,2)
            [~,dmrsAntIndices] = nrExtractResources(dmrsIndices(:,p),puschGrid);
            puschGrid(dmrsAntIndices) = puschGrid(dmrsAntIndices) + dmrsSymbols(:,p) * F(p,:);
        end

        % OFDM modulation
        txWaveform = hOFDMModulate(ue,puschGrid);

        % 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))]; %#ok<AGROW>
        [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 default behaviour of the channel model is to
        % apply this normalization to the received waveform
        SNR = 10^(SNRdB/20);
        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);
            [offset,mag] = nrPerfectTimingEstimate(pathGains,pathFilters);
        else
            % Practical synchronization. Correlate the received waveform
            % with the PUSCH 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(rxWaveform,ue.NRB,ue.SubcarrierSpacing,pusch.NSlot,dmrsIndices,dmrsSymbols,'CyclicPrefix',ue.CyclicPrefix); %#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 = hOFDMDemodulate(ue,rxWaveform);
        [K,L,R] = size(rxGrid);
        if (L < waveformInfo.SymbolsPerSlot)
            rxGrid = cat(2,rxGrid,zeros(K,waveformInfo.SymbolsPerSlot-L,R));
        end

        if (perfectChannelEstimator)
            % Perfect channel estimation, use the value of the path gains
            % provided by the channel
            estChannelGrid = nrPerfectChannelEstimate(pathGains,pathFilters,ue.NRB,ue.SubcarrierSpacing,pusch.NSlot,offset,sampleTimes,ue.CyclicPrefix);

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

            % Apply MIMO deprecoding to estChannelGrid to give an estimate
            % per transmission layer
            K = size(estChannelGrid,1);
            estChannelGrid = reshape(estChannelGrid,K*waveformInfo.SymbolsPerSlot*nRxAnts,nTxAnts);
            estChannelGrid = estChannelGrid * F.';
            if (strcmpi(pusch.TxScheme,'codebook'))
                W = nrPUSCHCodebook(pusch.NLayers,pusch.NAntennaPorts,pusch.TPMI,pusch.TransformPrecoding);
                estChannelGrid = estChannelGrid * W.';
            end
            estChannelGrid = reshape(estChannelGrid,K,waveformInfo.SymbolsPerSlot,nRxAnts,[]);
        else
            % Practical channel estimation between the received grid and
            % each transmission layer, using the PUSCH DM-RS for each layer
            [~,dmrsLayerIndices,dmrsLayerSymbols] = hPUSCHResources(ue,setfield(pusch,'TxScheme','nonCodebook')); %#ok<UNRCH>
            [estChannelGrid,noiseEst] = nrChannelEstimate(rxGrid,dmrsLayerIndices,dmrsLayerSymbols,'CyclicPrefix',ue.CyclicPrefix,'CDMLengths',puschIndicesInfo.CDMLengths);
        end

        % Get PUSCH resource elements from the received grid
        [puschRx,puschHest] = nrExtractResources(puschIndices,rxGrid,estChannelGrid);

        % Equalization
        [puschEq,csi] = nrEqualizeMMSE(puschRx,puschHest,noiseEst);

        % Decode PUSCH physical channel
        [ulschLLRs,rxSymbols] = nrPUSCHDecode(puschEq,pusch.Modulation,ue.NCellID,pusch.RNTI,noiseEst,pusch.TransformPrecoding,MRB);

        % Apply channel state information (CSI) produced by the equalizer,
        % including the effect of transform precoding if enabled
        if (pusch.TransformPrecoding)
            MSC = MRB * 12;
            csi = nrTransformDeprecode(csi,MRB) / sqrt(MSC);
            csi = repmat(csi((1:MSC:end).'),1,MSC).';
            csi = reshape(csi,size(rxSymbols));
        end
        csi = nrLayerDemap(csi);
        Qm = length(ulschLLRs) / length(rxSymbols);
        csi = reshape(repmat(csi{1}.',Qm,1),[],1);
        ulschLLRs = ulschLLRs .* csi;

        % Decode the UL-SCH transport channel
        decodeULSCH.TransportBlockLength = TBS;
        [decbits,harqProcesses(harqProcIdx).blkerr] = decodeULSCH(ulschLLRs,pusch.Modulation,pusch.NLayers,harqProcesses(harqProcIdx).RV,harqProcIdx-1);

        % Store values to calculate throughput
        simThroughput(snrIdx) = simThroughput(snrIdx) + (~harqProcesses(harqProcIdx).blkerr * TBS);
        maxThroughput(snrIdx) = maxThroughput(snrIdx) + TBS;

        % Display transport block error information
        if (displaySimulationInformation)
            fprintf('\n(%3.2f%%) HARQ Proc %d: ',100*(ue.NSymbol+size(puschGrid,2))/NSymbols,harqProcIdx);
            estrings = {'passed','failed'};
            rvi = harqProcesses(harqProcIdx).RVIdx;
            if rvi == 1
                ts = sprintf('Initial transmission (RV=%d)',rvSeq(rvi));
            else
                ts = sprintf('Retransmission #%d (RV=%d)',rvi-1,rvSeq(rvi));
            end
            fprintf('%s %s. ',ts,estrings{1+harqProcesses(harqProcIdx).blkerr});
        end

        % Update starting symbol number of next PUSCH transmission
        ue.NSymbol = ue.NSymbol + size(puschGrid,2);

        % Update count of overall number of PUSCH transmissions
        pusch.NSlot = pusch.NSlot + 1;

        % Update HARQ process counter
        harqProcCntr = harqProcCntr + 1;

    end

    % Display the results dynamically in the command window
    if (displaySimulationInformation)
        fprintf('\n');
    end
    fprintf([['\nThroughput(Mbps) for ' num2str(ue.NFrames) ' frame(s) '],'= %.4f\n'], 1e-6*simThroughput(snrIdx)/(ue.NFrames*10e-3));
    fprintf(['Throughput(%%) for ' num2str(ue.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;
if (pusch_init.TransformPrecoding)
    ofdmType = 'DFT-s-OFDM';
else
    ofdmType = 'CP-OFDM';
end
title(sprintf('%s / NRB=%d / SCS=%dkHz / %s %d/1024 / %dx%d', ...
    ofdmType,ue_init.NRB,ue_init.SubcarrierSpacing, ...
    pusch_init.Modulation, ...
    round(pusch_init.TargetCodeRate*1024),nTxAnts,nRxAnts));

simResults.simParameters = simParameters;
simResults.simThroughput = simThroughput;
simResults.maxThroughput = maxThroughput;