NFrames = 2;        % Number of frames

SNRIn = [-4 1];     % SNR range in dB

simulationParameters = []; % clear simulationParameters
simulationParameters.NDLRB = 50;
simulationParameters.PDSCH.TargetCodeRate = 0.5;
simulationParameters.PDSCH.Modulation = {'16QAM'};
simulationParameters.PDSCH.PRBSet = (0:9)';

txMode = 'TM8';     % TM7, TM8, TM9, TM10
DuplexMode = 'TDD'; % 'FDD', 'TDD'

switch txMode
    case 'TM7' % single layer (Port 5)
        simulationParameters.PDSCH.DCIFormat = 'Format1';
        simulationParameters.PDSCH.TxScheme = 'Port5';
        simulationParameters.PDSCH.NLayers = 1;
        ntxants = 4;
    case 'TM8' % up to two layers (dual layer Port 7-8, or single antenna
               % port, port 7 or 8)
        simulationParameters.PDSCH.DCIFormat = 'Format2B';
        simulationParameters.PDSCH.TxScheme = 'Port7-8'; %'Port7-8','Port8'
        simulationParameters.PDSCH.NLayers = 2;
        ntxants = 4;
        simulationParameters.PDSCH.NSCID = 0;
    case 'TM9' % up to eight layers (up to eight layers Port 7-14 or single
               % antenna port, port 7 or 8)
        simulationParameters.PDSCH.DCIFormat = 'Format2C';
        simulationParameters.PDSCH.TxScheme = 'Port7-14';
        simulationParameters.PDSCH.NLayers = 4;
        ntxants = 8;
        simulationParameters.CSIRefP = ntxants;
    case 'TM10' % up to eight layers (up to eight layers Port 7-14 or
                % single antenna port, port 7 or 8)
        simulationParameters.PDSCH.DCIFormat = 'Format2D';
        simulationParameters.PDSCH.TxScheme = 'Port7-14';
        simulationParameters.PDSCH.NLayers = 4;
        ntxants = 8;
        simulationParameters.CSIRefP = ntxants;
    otherwise
        error('Transmission mode should be one of TM7, TM8, TM9 or TM10.')
end

% For 2 codewords (more than 4 layers), the modulation field needs two values, one per codeword
if (simulationParameters.PDSCH.NLayers > 4) && (length(simulationParameters.PDSCH.Modulation)<2)
    simulationParameters.PDSCH.Modulation = repmat(simulationParameters.PDSCH.Modulation,1,2);
end

% Initialize W to zero
simulationParameters.PDSCH.W = zeros(simulationParameters.PDSCH.NLayers, ntxants);

simulationParameters.DuplexMode = DuplexMode;
simulationParameters.TotSubframes = 1;

enb = lteRMCDL(simulationParameters);

trBlkSizes = enb.PDSCH.TrBlkSizes;
rvSequence = enb.PDSCH.RVSeq;

ncw = length(string(enb.PDSCH.Modulation));

hDisplayENBParameterSummary(enb, txMode);

channel = struct;                    % Channel config structure
channel.Seed = 6;                    % Channel seed
channel.NRxAnts = 8;                 % Number of receive antennas
channel.DelayProfile ='EPA';         % Delay profile
channel.DopplerFreq = 5;             % Doppler frequency
channel.MIMOCorrelation = 'Low';     % Multi-antenna correlation
channel.NTerms = 16;                 % Oscillators used in fading model
channel.ModelType = 'GMEDS';         % Rayleigh fading model type
channel.InitPhase = 'Random';        % Random initial phases
channel.NormalizePathGains = 'On';   % Normalize delay profile power
channel.NormalizeTxAnts = 'On';      % Normalize for transmit antennas
% The channel sampling rate depends on the FFT size used in the OFDM
% modulator. This can be obtained using the function lteOFDMInfo.
ofdmInfo = lteOFDMInfo(enb);
channel.SamplingRate = ofdmInfo.SamplingRate;

cec = struct;                        % Channel estimation config structure
cec.PilotAverage = 'UserDefined';    % Type of pilot symbol averaging
cec.FreqWindow = 41;                 % Frequency window size
cec.TimeWindow = 27;                 % Time window size
cec.InterpType = 'Cubic';            % 2D interpolation type
cec.InterpWindow = 'Centered';       % Interpolation window type
cec.InterpWinSize = 1;               % Interpolation window size
cec.Reference = 'DMRS';              % Use DMRS as reference signal

displaySimulationInformation = true;

% Initialize variables used in the simulation and analysis
% 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);
% Array to store RV sequence index history
allRvSeqIdxHistory = cell(1,numel(SNRIn));
% Copy of channel structure to optimize parallel processing (only if
% running the example with Parallel Computing Toolbox)
channelInit = channel;
% During simulation some fields of enb will be updated, make a copy to
% reinitialize it when simulating each SNR point
enbInit = enb;

% For TDD precalculate vector of subframe types: D, S and U for downlink,
% special and uplink respectively
if strcmpi(DuplexMode,'TDD')
    subframeType = char(10,1);
    initialSubframeNo = enb.NSubframe;
    for sNo=0:9 % for all subframes in a frame
        enb.NSubframe = sNo;
        duplexInfo = lteDuplexingInfo(enb);
        subframeType(sNo+1) = duplexInfo.SubframeType(1); % first char: D, S or U
    end
    enb.NSubframe = initialSubframeNo;
end

% CFI can be a scalar or a vector, create a local copy of CFI as a vector
% (one value per subframe)
if numel(enb.CFI) == 1
    CFI = repmat(enb.CFI,1,10);
else
    CFI = enb.CFI;
end

% Get the HARQ ID sequence for HARQ processing. This is a list of indices
% for HARQ process scheduling.
[~,~,enbOut] = lteRMCDLTool(enb, []);
harqProcessSequenceInit = enbOut.PDSCH.HARQProcessSequence;

for snrIdx = 1:numel(SNRIn)     % comment out for parallel computing
%parfor snrIdx = 1:numel(SNRIn) % uncomment for parallel computing

    % Set the random number generator seed depending on the loop variable
    % to ensure independent random streams
    rng(snrIdx,'combRecursive');

    % Reinitialize variables (they are modified during each SNR point
    % simulation)
    enb = enbInit;
    channel = channelInit;
    harqProcessSequence = harqProcessSequenceInit;

    % Initialize the state of all HARQ processes
    harqProcesses = hNewHARQProcess(enb);

    % Set up variables for the main loop
    lastOffset = 7;       % Initialize overall frame timing offset
                          % (set equal to channel implementation delay)
    frameOffset = 7;      % Initialize frame timing offset
                          % (set equal to channel implementation delay)
    blkCRC = [];          % Block CRC for all considered subframes
    bitTput = [];         % Number of successfully received bits per subframe
    txedTrBlkSizes = [];  % Number of transmitted bits per subframe
    precodingMatrix = []; % Precoding matrix
    rxSymbols = [];       % DL-SCH symbols for constellation plot (codeword 1)
    rxSymbols2 = [];      % DL-SCH symbols for constellation plot (codeword 2)

    % The variable rvSeqIdxHistory stores the history of the value of the
    % RV index sequence for all the HARQ processes
    rvSeqIdxHistory = nan(ncw, NFrames*10);

    % Flag to indicate if a precoding matrix W is available. This is useful
    % if no W has been calculated at the beginning of the simulation.
    isWready = false;
    % Flag to indicate is a subframe is to be processed. Set to true if
    % there is data to be processed in the subframe, i.e. DL subframe or
    % non-zero transport block size.
    processSubframe = false;

    % Main for loop: for all subframes
    for subframeNo = 0:(NFrames*10-1)

        % Update subframe number
        enb.NSubframe = subframeNo;

        % Load CFI for current subframe
        enb.CFI = CFI(mod(subframeNo,length(CFI))+1);

        % Channel time for the current subframe
        channel.InitTime = subframeNo/1000;

        % Get HARQ process ID for the subframe from HARQ process sequence
        harqID = harqProcessSequence(mod(subframeNo, length(harqProcessSequence))+1);

        % Extract the current subframe transport block size(s)
        trBlk = trBlkSizes(:, mod(subframeNo, 10)+1).';

        % Load precoding matrix calculated in previous subframe if it
        % exists and if there is data to transmit. Set the flag to trigger
        % subframe processing.
        if isWready && any(trBlk)
            harqProcesses(harqID).txConfig.W = precodingMatrix;
            processSubframe = true;
        else
            processSubframe = false;
        end

        % Precoding matrix calculation
        if strcmpi(DuplexMode,'TDD')
            % Estimate channel in UL subframe
            if strcmp(subframeType(mod(subframeNo,10)+1),'U')
                processSubframe = false; % UL subframe, no DL data
                % Only perform channel estimate if next subframe is DL
                if strcmp(subframeType(mod((subframeNo+1),10)+1),'D')
                    precodingMatrix = hCalculatePrecodingMatrix(enb, channel);
                    isWready = true;
                end
            end
        else %FDD
            % Get transport block for next subframe
            trBlkNext = trBlkSizes(:, mod(subframeNo+1, 10)+1).';

            % Calculate the precoding matrix for next subframe only if it
            % carries data (i.e. non-zero trBlkNext)
            if any(trBlkNext)
                precodingMatrix = hCalculatePrecodingMatrix(enb, channel);
                isWready = true;
            else
                isWready = false;
            end
        end

        % Subframe processing
        if processSubframe

            % Update current HARQ process
            harqProcesses(harqID) = hHARQScheduling( ...
                harqProcesses(harqID), subframeNo, rvSequence);

            % Display run time information
            if displaySimulationInformation && any(trBlk)
                disp(' ');
                disp(['Subframe: ' num2str(subframeNo) '. HARQ process index: ' num2str(harqID)]);
            end

            % Update RV sequence index table
            rvSeqIdxHistory(:,subframeNo+1) = ...
                harqProcesses(harqID).txConfig.RVIdx.';

            % Update the PDSCH transmission configuration with HARQ
            % process state
            enb.PDSCH = harqProcesses(harqID).txConfig;

            % Cell payload
            dlschTransportBlk = harqProcesses(harqID).data;

            % Create transmit waveform
            txWaveform = lteRMCDLTool(enb, dlschTransportBlk);

            % Add 25 sample padding. This is to cover the range of delays
            % expected from channel modeling (a combination of
            % implementation delay and channel delay spread)
            txWaveform =  [txWaveform; zeros(25, ntxants)]; %#ok<AGROW>

            % Pass data through channel model
            rxNoiselessWaveform = lteFadingChannel(channel,txWaveform);

            % Calculate noise gain including compensation for downlink
            % power allocation
            SNR = 10^((SNRIn(snrIdx)-enb.PDSCH.Rho)/20);

            % Normalize noise power to take account of sampling rate,
            % which is a function of the IFFT size used in OFDM
            % modulation, and the number of antennas
            N0 = 1/(sqrt(2.0*ntxants*double(ofdmInfo.Nfft))*SNR);

            % Create additive white Gaussian noise
            noise = N0*complex(randn(size(rxNoiselessWaveform)),...
                randn(size(rxNoiselessWaveform)));

            % Add AWGN to the received time domain waveform and scale
            % for required power
            rxWaveform = rxNoiselessWaveform + noise;

            % Receiver
            % Once every frame, on subframe 0, calculate a new
            % synchronization offset. An offset within the range of
            % delays expected from the channel modeling (a combination
            % of implementation delay and channel delay spread)
            % indicates success
            if (mod(subframeNo,10)==0)
                frameOffset = lteDLFrameOffset(enb,rxWaveform);
                if (frameOffset > 25)
                    frameOffset = lastOffset;
                end
                lastOffset = frameOffset;
            end
            rxWaveform = rxWaveform(1+frameOffset:end,:);

            % Perform OFDM demodulation on the received data to obtain
            % the resource grid
            rxSubframe = lteOFDMDemodulate(enb,rxWaveform);

            % Channel estimation
            [estChannelGrid,noiseEst] = ...
                lteDLChannelEstimate(enb,enb.PDSCH,cec,rxSubframe);

            % Perform equalization, deprecoding, layer demapping,
            % demodulation and descrambling on the received data using the
            % channel estimate.

            % Get PDSCH indices
            pdschIndices = ltePDSCHIndices(enb,enb.PDSCH,enb.PDSCH.PRBSet);

            % Get PDSCH resource elements. Scale the received subframe
            % by the PDSCH power factor Rho. The PDSCH is scaled by
            % Rho, while the cell reference symbols used for channel
            % estimation (used in the PDSCH decoding stage) are not.
            [pdschRx, pdschHest] = lteExtractResources(pdschIndices, ...
                rxSubframe*(10^(-enb.PDSCH.Rho/20)),estChannelGrid);

            % Decode PDSCH
            [dlschBits,dlschSymbols] = ltePDSCHDecode(enb,enb.PDSCH,...
                pdschRx,pdschHest,noiseEst);

            % Store the decoded DLSCH symbols for constellation
            % plotting
            rxSymbols = [rxSymbols; dlschSymbols{1}(:)]; %#ok<AGROW>
            if ncw>1
                rxSymbols2 = [rxSymbols2; dlschSymbols{2}(:)]; %#ok<AGROW>
            end

            % Decode the DL-SCH
            [decbits,harqProcesses(harqID).blkerr,harqProcesses(harqID).decState] = ...
                lteDLSCHDecode(enb,enb.PDSCH,trBlk,dlschBits, ...
                harqProcesses(harqID).decState);

            % Display block errors
            if displaySimulationInformation
                if any(harqProcesses(harqID).blkerr)
                    disp(['Block error. RV index: ' num2str(harqProcesses(harqID).txConfig.RVIdx) ', CRC: ' num2str(harqProcesses(harqID).blkerr)])
                else
                    disp(['No error. RV index: ' num2str(harqProcesses(harqID).txConfig.RVIdx) ', CRC: ' num2str(harqProcesses(harqID).blkerr)])
                end
            end

            % Store values needed to calculate throughput
            % Only for subframes with data
            if(any(trBlk))
                blkCRC = [blkCRC harqProcesses(harqID).blkerr]; %#ok<AGROW>
                bitTput = [bitTput trBlk.*(1- ...
                    harqProcesses(harqID).blkerr)]; %#ok<AGROW>
                txedTrBlkSizes = [txedTrBlkSizes trBlk]; %#ok<AGROW>
            end
        end
    end

    % Plot received symbol constellation
    if displaySimulationInformation
        figure;plot(rxSymbols,'.r'); title(['Received constellation. SNR ' ...
                                 num2str(SNRIn(snrIdx)) ' dB. Codeword 1'])
        xlabel('In-phase Amplitude');
        ylabel('Quadrature Amplitude');
        if ncw>1
            figure;plot(rxSymbols2,'.r'); title(['Received constellation. SNR ' ...
                                 num2str(SNRIn(snrIdx)) ' dB. Codeword 2'])
            xlabel('In-phase Amplitude');
            ylabel('Quadrature Amplitude');
        end
    end

    % Calculate the 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('\nSNR = %.2f dB. Throughput for %d Frame(s) = %.4f Mbps\n',...
         SNRIn(snrIdx),NFrames,1e-6*simThroughput(snrIdx)/(NFrames*10e-3));
    fprintf('SNR = %.2f dB. Throughput(%%) for %d Frame(s) = %.4f %%\n',...
    SNRIn(snrIdx),NFrames,simThroughput(snrIdx)*100/maxThroughput(snrIdx));

    allRvSeqIdxHistory{snrIdx} = rvSeqIdxHistory;

end

hPlotRVSequence(SNRIn,allRvSeqIdxHistory,NFrames);

% Plot throughput
figure
legendString = [char(txMode) ': ' num2str(simulationParameters.PDSCH.NLayers) ...
                ' layer(s), ' num2str(ntxants) ' ant(s)'];
plot(SNRIn, simThroughput*100./maxThroughput,'*-.');
xlabel('SNR (dB)');
ylabel('Throughput (%)');
title('Throughput vs SNR');
legend(legendString,'Location','NorthWest');
grid on;