NFrames = 2;                % Number of frames
SNRIn = [10.3 12.3 14.3];   % SNR range in dB

txMode = 'TM4'; % TM1, TM2, TM3, TM4, TM6

simulationParameters = []; % clear simulationParameters
simulationParameters.NDLRB = 50;
simulationParameters.PDSCH.TargetCodeRate = 0.5;
simulationParameters.PDSCH.PRBSet = (0:49)';

switch txMode
% Single antenna (Port0) mode (TM1)
    case 'TM1'
        fprintf('\nTM1 - Single antenna (port 0)\n');
        simulationParameters.PDSCH.TxScheme = 'Port0';
        simulationParameters.PDSCH.DCIFormat = 'Format1';
        simulationParameters.CellRefP = 1;
        simulationParameters.PDSCH.Modulation = {'16QAM'};

% Transmit diversity mode (TM2)
    case 'TM2'
        fprintf('\nTM2 - Transmit diversity\n');
        simulationParameters.PDSCH.TxScheme = 'TxDiversity';
        simulationParameters.PDSCH.DCIFormat = 'Format1';
        simulationParameters.CellRefP = 2;
        simulationParameters.PDSCH.Modulation = {'16QAM'};
        simulationParameters.PDSCH.NLayers = 2;

% CDD mode (TM3)
    case 'TM3'
        fprintf('\nTM3 - CDD\n');
        simulationParameters.PDSCH.TxScheme = 'CDD';
        simulationParameters.PDSCH.DCIFormat = 'Format2A';
        simulationParameters.CellRefP = 2;
        simulationParameters.PDSCH.Modulation = {'16QAM', '16QAM'};
        simulationParameters.PDSCH.NLayers = 2;

% Spatial multiplexing mode (TM4)
    case 'TM4'
        fprintf('\nTM4 - Codebook based spatial multiplexing\n');
        simulationParameters.CellRefP = 2;
        simulationParameters.PDSCH.Modulation = {'16QAM', '16QAM'};
        simulationParameters.PDSCH.DCIFormat = 'Format2';
        simulationParameters.PDSCH.TxScheme = 'SpatialMux';
        simulationParameters.PDSCH.NLayers = 2;
        % No codebook restriction
        simulationParameters.PDSCH.CodebookSubset = '';

% Single layer spatial multiplexing mode (TM6)
    case 'TM6'
        fprintf(...
        '\nTM6 - Codebook based spatial multiplexing with single layer\n');
        simulationParameters.CellRefP = 4;
        simulationParameters.PDSCH.Modulation = {'QPSK'};
        simulationParameters.PDSCH.DCIFormat = 'Format2';
        simulationParameters.PDSCH.TxScheme = 'SpatialMux';
        simulationParameters.PDSCH.NLayers = 1;
        % No codebook restriction
        simulationParameters.PDSCH.CodebookSubset = '';

    otherwise
        error('Transmission mode should be one of TM1, TM2, TM3, TM4 or TM6.')
end

% Set other simulationParameters fields applying to all TMs
simulationParameters.TotSubframes = 1; % Generate one subframe at a time
simulationParameters.PDSCH.CSI = 'On'; % Soft bits are weighted by CSI

enb = lteRMCDL(simulationParameters);

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

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

pmiDelay = 8;

hDisplayENBParameterSummary(enb, txMode);

channel.Seed = 6;                    % Channel seed
channel.NRxAnts = 2;                 % 2 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 sampling rate for the channel model is set using the value returned
% from <matlab:doc('lteOFDMInfo') lteOFDMInfo>.

ofdmInfo = lteOFDMInfo(enb);
channel.SamplingRate = ofdmInfo.SamplingRate;

% Perfect channel estimator flag
perfectChanEstimator = false;

% Configure channel estimator
cec.PilotAverage = 'UserDefined';   % Type of pilot symbol averaging
cec.FreqWindow = 41;                % Frequency window size in REs
cec.TimeWindow = 27;                % Time window size in REs
cec.InterpType = 'Cubic';           % 2D interpolation type
cec.InterpWindow = 'Centered';      % Interpolation window type
cec.InterpWinSize = 1;              % Interpolation window size

displaySimulationInformation = true;

% The number of transmit antennas P is obtained from the resource grid
% dimensions. 'dims' is M-by-N-by-P where M is the number of subcarriers, N
% is the number of symbols and P is the number of transmit antennas.
dims = lteDLResourceGridSize(enb);
P = dims(3);

% 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);

% The temporary variables 'enb_init' and 'channel_init' are used to create
% the temporary variables 'enb' and 'channel' within the SNR loop to create
% independent simulation loops for the 'parfor' loop
enb_init = enb;
channel_init = channel;
legendString = ['Throughput: ' char(enb.PDSCH.TxScheme)];
allRvSeqPtrHistory = cell(1,numel(SNRIn));
nFFT = ofdmInfo.Nfft;

for snrIdx = 1:numel(SNRIn)
% parfor snrIdx = 1:numel(SNRIn)
% To enable the use of parallel computing for increased speed comment out
% the 'for' statement above and uncomment the 'parfor' statement below.
% This needs the Parallel Computing Toolbox. If this is not installed
% 'parfor' will default to the normal 'for' statement. If 'parfor' is
% used it is recommended that the variable 'displaySimulationInformation'
% above is set to false, otherwise the simulation information displays for
% each SNR point will overlap.

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

    SNRdB = SNRIn(snrIdx);
    fprintf('\nSimulating at %g dB SNR for %d Frame(s)\n' ,SNRdB, NFrames);

    % Initialize variables used in the simulation and analysis
    offsets = 0;            % Initialize frame offset value
    offset = 0;             % Initialize frame offset value for radio frame
    blkCRC = [];            % Block CRC for all considered subframes
    bitTput = [];           % Number of successfully received bits per subframe
    txedTrBlkSizes = [];    % Number of transmitted bits per subframe
    enb = enb_init;         % Initialize RMC configuration
    channel = channel_init; % Initialize channel configuration
    pmiIdx = 0;             % PMI index in delay queue

    % The variable harqPtrTable stores the history of the value of the
    % pointer to the RV sequence values for all the HARQ processes.
    % Pre-allocate with NaNs as some subframes do not have data
    rvSeqPtrHistory = NaN(ncw, NFrames*10);

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

    % Use random PMIs for the first 'pmiDelay' subframes until feedback is
    % available from the UE; note that PMI feedback is only applicable for
    % spatial multiplexing TMs (TM4 and TM6), but the code here is required
    % for complete initialization of variables in the SNR loop when using
    % the Parallel Computing Toolbox.
    pmidims = ltePMIInfo(enb,enb.PDSCH);
    txPMIs = randi([0 pmidims.MaxPMI], pmidims.NSubbands, pmiDelay);
    % Initialize HARQ process IDs to 1 as the first non-zero transport
    % block will always be transmitted using the first HARQ process. This
    % will be updated with the full sequence output by lteRMCDLTool after
    % the first call to the function
    harqProcessSequence = 1;

    for subframeNo = 0:(NFrames*10-1)

        % Update subframe number
        enb.NSubframe = subframeNo;

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

        % If there is a transport block scheduled in the current subframe
        % (indicated by non-zero 'harqID'), perform transmission and
        % reception. Otherwise continue to the next subframe
        if harqID == 0
            continue;
        end

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

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

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

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

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

        % Set the PMI to the appropriate value in the delay queue
        if strcmpi(enb.PDSCH.TxScheme,'SpatialMux')
            pmiIdx = mod(subframeNo, pmiDelay);  % PMI index in delay queue
            enb.PDSCH.PMISet = txPMIs(:, pmiIdx+1); % Set PMI
        end

        % Create transmit waveform and get the HARQ scheduling ID sequence
        % from 'enbOut' structure output which also contains the waveform
        % configuration and OFDM modulation parameters
        [txWaveform,~,enbOut] = lteRMCDLTool(enb, data);

        % 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, P)]; %#ok<AGROW>

        % Get the HARQ ID sequence from 'enbOut' for HARQ processing
        harqProcessSequence = enbOut.PDSCH.HARQProcessSequence;

        % Initialize channel time for each subframe
        channel.InitTime = subframeNo/1000;

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

        % Calculate noise gain including compensation for downlink power
        % allocation
        SNR = 10^((SNRdB-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*enb.CellRefP*double(nFFT))*SNR);

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

        % Add AWGN to the received time domain waveform
        rxWaveform = rxWaveform + noise;

        % Once every frame, on subframe 0, calculate a new synchronization
        % offset
        if (mod(subframeNo,10) == 0)
            offset = lteDLFrameOffset(enb, rxWaveform);
            if (offset > 25)
                offset = offsets(end);
            end
            offsets = [offsets offset]; %#ok
        end

        % Synchronize the received waveform
        rxWaveform = rxWaveform(1+offset:end, :);

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

        % Channel estimation
        if(perfectChanEstimator)
            estChannelGrid = lteDLPerfectChannelEstimate(enb, channel, offset); %#ok
            noiseGrid = lteOFDMDemodulate(enb, noise(1+offset:end ,:));
            noiseEst = var(noiseGrid(:));
        else
            [estChannelGrid, noiseEst] = lteDLChannelEstimate( ...
                enb, enb.PDSCH, cec, rxSubframe);
        end

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

        % Get PDSCH resource elements from the received subframe. Scale the
        % received subframe by the PDSCH power factor Rho. The PDSCH is
        % scaled by this amount, 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 = ltePDSCHDecode(...
                             enb, enb.PDSCH, pdschRx, pdschHest, noiseEst);

        % 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 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

        % Provide PMI feedback to the eNodeB
        if strcmpi(enb.PDSCH.TxScheme,'SpatialMux')
            PMI = ltePMISelect(enb, enb.PDSCH, estChannelGrid, noiseEst);
            txPMIs(:, pmiIdx+1) = PMI;
        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([['\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));

    allRvSeqPtrHistory{snrIdx} = rvSeqPtrHistory;

end

% Plot the RV sequence for all HARQ processes
hPlotRVSequence(SNRIn,allRvSeqPtrHistory,NFrames);

% Plot throughput
figure
plot(SNRIn, simThroughput*100./maxThroughput,'*-.');
xlabel('SNR (dB)');
ylabel('Throughput (%)');
title('Throughput vs. SNR')
legend(legendString,'Location','NorthWest');
grid on;