NFrames = 1;                % Number of frames to simulate at each SNR
SNRIn = [-4.1, -2.0, 0.1];  % SNR points to simulate

ue.TotSubframes = 1;    % Total number of subframes to generate a waveform for
ue.NCellID = 10;        % Cell identity
ue.RC = 'A3-2';         % FRC number
ue.PUSCH.NLayers = 1;   % Number of layers
ue.NTxAnts = 1;         % Number of transmit antennas

chcfg.NRxAnts = 2;               % Number of receive antennas
chcfg.DelayProfile = 'EPA';      % Delay profile
chcfg.DopplerFreq = 5.0;         % Doppler frequency
chcfg.MIMOCorrelation = 'Low';   % MIMO correlation
chcfg.Seed = 100;                % Channel seed
chcfg.NTerms = 16;               % Oscillators used in fading model
chcfg.ModelType = 'GMEDS';       % Rayleigh fading model type
chcfg.InitPhase = 'Random';      % Random initial phases
chcfg.NormalizePathGains = 'On'; % Normalize delay profile power
chcfg.NormalizeTxAnts = 'On';    % Normalize for transmit antennas

if ue.PUSCH.NLayers == 1
   freqWindow = 13;
else
   freqWindow = 12; % Averaging to take advantage of orthogonal DMRS for multilayer
end
cec.FreqWindow = freqWindow;      % Frequency averaging windows in REs
cec.PilotAverage = 'UserDefined'; % Type of pilot averaging
cec.TimeWindow = 1;               % Time averaging windows in REs
cec.InterpType = 'cubic';         % Interpolation type
cec.Reference = 'Antennas';       % Reference for channel estimation

% Generate FRC configuration structure for A3-2
frc = lteRMCUL(ue);

% Transport block sizes for each subframe within a frame
trBlkSizes = frc.PUSCH.TrBlkSizes;
codedTrBlkSizes = frc.PUSCH.CodedTrBlkSizes;

% Redundancy version sequence
rvSequence = frc.PUSCH.RVSeq;

info = lteSCFDMAInfo(frc);
chcfg.SamplingRate = info.SamplingRate;

% Initialize variables used in the simulation and analysis
totalBLKCRC = zeros(numel(SNRIn), NFrames*10);   % Total block CRC vector
bitThroughput = zeros(numel(SNRIn), NFrames*10); % Total throughput vector
resultIndex = 1;        % Initialize frame counter index

for SNRdB = SNRIn

    fprintf('\nSimulating at %g dB SNR for a total %d Frame(s)', ...
        SNRdB, NFrames);

    % Calculate required AWGN channel noise
    SNR = 10^(SNRdB/20);
    N = 1/(SNR*sqrt(double(info.Nfft)))/sqrt(2.0);
    rng('default');

    % Store results for every subframe at SNR point
    bitTp = zeros(1, NFrames*10);  % Intermediate bit throughput vector
    blkCRC = zeros(1, NFrames*10); % Intermediate block CRC vector

    % Initialize state of all HARQ processes
    harqProcesses = hNewHARQProcess(frc);
    % 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 lteRMCULTool after
    % the first call to the function
    harqProcessSequence = 1;

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

        % Update subframe number
        frc.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);

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

        % Create transmit waveform and get the HARQ scheduling ID sequence
        % from 'frcOut' structure output which also contains the waveform
        % configuration and OFDM modulation parameters
        [txWaveform,~,frcOut] = lteRMCULTool(frc, 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, frc.NTxAnts)]; %#ok<AGROW>

        % Get the HARQ ID sequence from 'frcOut' for HARQ processing
        harqProcessSequence = frcOut.PUSCH.HARQProcessSequence;

        % The initialization time for channel modeling is set each subframe
        % to simulate a continuously varying channel
        chcfg.InitTime = subframeNo/1000;

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

        % Add noise at the receiver
        v = N*complex(randn(size(rxWaveform)), randn(size(rxWaveform)));
        rxWaveform = rxWaveform+v;

        % Calculate synchronization offset
        offset = lteULFrameOffset(frc, frc.PUSCH, rxWaveform);
        if (offset < 25)
            offsetused = offset;
        end

        % SC-FDMA demodulation
        rxSubframe = lteSCFDMADemodulate(frc, ...
            rxWaveform(1+offsetused:end, :));

        % Channel and noise power spectral density estimation
        [estChannelGrid, noiseEst] = lteULChannelEstimate(frc, ...
            frc.PUSCH, cec, rxSubframe);

        % Extract REs corresponding to the PUSCH from the given subframe
        % across all receive antennas and channel estimates
        puschIndices = ltePUSCHIndices(frc, frc.PUSCH);
        [puschRx, puschEstCh] = lteExtractResources( ...
            puschIndices, rxSubframe, estChannelGrid);

        % MMSE equalization and symbol scaling
        rxSymbols = lteEqualizeMMSE(puschRx, puschEstCh, noiseEst);
        rxSymbols = rxSymbols*(double(frc.NTxAnts/frc.PUSCH.NLayers));

        % Update frc.PUSCH to carry complete information of the UL-SCH
        % coding configuration
        trBlkSize = trBlkSizes(mod(subframeNo, 10)+1);
        frc.PUSCH = lteULSCHInfo(frc, frc.PUSCH, trBlkSize, 'chsconcat');

        % Decode the PUSCH
        rxEncodedBits = ltePUSCHDecode(frc, frc.PUSCH, rxSymbols);

        % Decode the UL-SCH channel and store the block CRC error for given
        % HARQ process
        [rxDecodedBits, harqProcesses(harqID).blkerr, ...
            harqProcesses(harqID).decState] = lteULSCHDecode(...
            frc, frc.PUSCH, trBlkSize, ...
            rxEncodedBits, harqProcesses(harqID).decState);

        % Store the CRC calculation and total number of bits per subframe
        % successfully decoded
        blkCRC(subframeNo+1) = harqProcesses(harqID).blkerr;
        bitTp(subframeNo+1) = trBlkSize.*(1-harqProcesses(harqID).blkerr);

    end

    % Record the block CRC error and bit throughput for the total number of
    % frames simulated at a particular SNR
    totalBLKCRC(resultIndex, :) = blkCRC;
    bitThroughput(resultIndex, :) = bitTp;
    resultIndex = resultIndex + 1;
end

% Throughput calculation as a percentage
throughput = 100*(1-mean(totalBLKCRC, 2)).';

hPUSCHResults(SNRIn, NFrames, trBlkSizes, throughput, bitThroughput);