NFrames = 10; % Number of frames
SNRIn = -0.1; % SNR range in dB

frc.RC = 'A3-2';                 % FRC number
frc.DuplexMode = 'FDD';          % Duplex mode
frc.TotSubframes = 1;            % Total number of subframes to generate
frc.NCellID = 10;                % Cell identity
frc.CyclicPrefixUL = 'Normal';   % Uplink cyclic prefix length
frc.CyclicPrefix = 'Normal';     % Downlink cyclic prefix length
frc.NTxAnts = 4;                 % Number of transmit antennas

frc.SRS.SubframeConfig = 1;    % Cell-specific schedule: 2 ms periodicity
frc.SRS.NTxAnts = frc.NTxAnts; % Configure SORTD same as no of UE antennas
frc.SRS.BWConfig = 7;          % Cell-specific SRS bandwidth
frc.SRS.BW = 0;                % UE-specific SRS bandwidth
frc.SRS.ConfigIdx = 0;         % UE-specific schedule:  2ms periodicity
frc.SRS.CyclicShift = 0;       % Cyclic shift 0
frc.SRS.SeqGroup = 0;          % Sequence group 0
frc.SRS.SeqIdx = 0;            % Base sequence 0
frc.SRS.TxComb = 0;            % Transmission comb 0
frc.SRS.FreqPosition = 0;      % Frequency-domain position 0
frc.SRS.HoppingBW = 0;         % Disable hopping (as HoppingBW >=0 BW)

pmiDelay = 8;

chcfg.NRxAnts = 4;               % Number of receive antennas
chcfg.DelayProfile = 'EPA';      % Delay profile
chcfg.DopplerFreq = 5.0;         % Doppler frequency
chcfg.MIMOCorrelation = 'Low';   % MIMO correlation
chcfg.Seed = 9;                  % 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

% Controls channel estimator behavior
perfectChanEstimator = false;      % Valid values are true or false

% Channel estimator configuration
cec.PilotAverage = 'UserDefined'; % Type of pilot symbol averaging
cec.FreqWindow = 12;              % Frequency averaging window in REs (special mode)
cec.TimeWindow = 1;               % Time averaging window in REs (special mode)
cec.InterpType = 'Cubic';         % 2D interpolation type

% Point of reference for the imperfect channel estimator
reference = 'Antennas';

% Validate reference is present when imperfect channel estimation used
if (~perfectChanEstimator)
    if (~any(strcmpi(reference,{'Antennas' 'Layers'})))
        error(['Channel estimator reference must be either ''Antennas'...
        ' or ''Layers''.']);
    end
else
    reference = 'Antennas';
end

% Generate parameter structure for single codeword A3-2 FRC
frc = lteRMCUL(frc);

% Then update Physical Uplink Shared Channel (PUSCH) parameters to define
% two identically configured single layer codewords
frc.PUSCH.NLayers = 2;  % the layers are shared across the 2 codewords
frc.PUSCH.Modulation = repmat({frc.PUSCH.Modulation}, 1, 2);
frc.PUSCH.TrBlkSizes = repmat(frc.PUSCH.TrBlkSizes, 2, 1);
frc.PUSCH.RVSeq = repmat(frc.PUSCH.RVSeq, 2, 1);

% Record transport block sizes for each subframe in a frame
trBlkSizes = frc.PUSCH.TrBlkSizes;

% The number of codewords is the number of transport block sizes
ncw = size(trBlkSizes,1);

% Record RV sequence
rvSequence = frc.PUSCH.RVSeq;

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

% Store results for each SNR point and each subframe containing data for
% the whole simulation
nDataTBS = sum(trBlkSizes(:)~=0)*NFrames;
crc = zeros(numel(SNRIn), nDataTBS);  % Total block CRC error vector
tput = zeros(numel(SNRIn), nDataTBS); % Total throughput vector
ResultIndex = 1; % SNR point result index

for SNRdB = SNRIn

    % Configure random number generators
    rng('default');

    % Noise configuration
    fprintf('\nSimulating at %gdB SNR for a total %d Frame(s)', ...
            SNRdB, NFrames);
    SNR = 10^(SNRdB/20);
    N = 1/(SNR*sqrt(double(dims.Nfft)))/sqrt(2.0*frc.NTxAnts);

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

    % Use random PMIs for the first 'pmiDelay' subframes until feedback is
    % available from the UE if multiple transmit antennas are used
    if (frc.NTxAnts>1)
        pmidims = lteULPMIInfo(frc, frc.PUSCH);
        txPMIs = randi([0 pmidims.MaxPMI], pmidims.NSubbands, pmiDelay);
    end

    % Initialize result store for SNR point tested
    crcSNR = zeros(nDataTBS/ncw, ncw);  % Intermediate block CRC
    tputSNR = zeros(nDataTBS/ncw, ncw); % Intermediate throughput
    dataSubframeIndex = 1;

    % Loop for all subframes at this SNR
    offsetused = 0;
    for subframeNo = 0:(NFrames*10-1)

        % Update subframe number
        frc.NSubframe = mod(subframeNo, 10);

        % If this is an uplink subframe
        duplexDims = lteDuplexingInfo(frc);
        if(strcmp(duplexDims.SubframeType,'Uplink')==1)

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

            % Set the PMI to the appropriate value in the delay queue (if
            % multiple antennas are being used)
            if (frc.NTxAnts>1)
                pmiIdx = mod(subframeNo, pmiDelay);
                frc.PUSCH.PMI = txPMIs(:, pmiIdx+1);
            end

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

            % Pass data through the fading channel model. The
            % initialization time for channel modeling is set each subframe
            % to simulate a continuously varying channel.
            chcfg.InitTime = subframeNo/1000;
            rxWaveform = lteFadingChannel(chcfg, txWaveform);

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

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

            % Mark this subframe as not being involved in PMI update
            updatePMI = false;

            % Channel and noise estimation
            if (perfectChanEstimator)
                % Create reference grid - not required for perfect channel
                % estimator, but needed as an input to lteULPMISelect to
                % match call with imperfect channel estimator.
                refGrid = NaN(lteULResourceGridSize(frc)); %#ok

                % Perfect channel estimation, perfect knowledge of all REs
                % and perfect knowledge of channel noise
                estChannelGrid = lteULPerfectChannelEstimate( ...
                    frc, chcfg, offsetused);
                n = lteSCFDMADemodulate(frc, noise(1+offsetused:end, :));
                noiseEst = var(reshape(n, numel(n), 1));
                updatePMI = true;
            else
                % Use imperfect channel estimator with support from the SRS
                % if available
                if (strcmpi(reference, 'Antennas')==1 || ...
                        (frc.NTxAnts == frc.PUSCH.NLayers))
                    refGrid = NaN(lteULResourceGridSize(frc));

                    if (isfield(frc, 'SRS'))
                        [srsIndices, srsIndicesDims] = lteSRSIndices( ...
                            frc, frc.SRS);
                        if (~isempty(srsIndices))
                            srsSymbols = lteSRS(frc, frc.SRS);
                            refGrid(srsIndices) = srsSymbols;
                            updatePMI = true;
                        end
                    end
                elseif (strcmpi(reference, 'Layers')==1)
                    refGrid = NaN(lteULResourceGridSize ...
                              (frc, frc.PUSCH.NLayers));
                end
                cec.Reference = reference;
                [estChannelGrid, noiseEst] = lteULChannelEstimate( ...
                    frc, frc.PUSCH, cec, rxSubframe, refGrid);
            end

            % Shorten PUSCH capacity shortening as appropriate
            if (isfield(frc, 'SRS'))
                srsDims = lteSRSInfo(frc, frc.SRS);
                frc.Shortened = srsDims.IsSRSSubframe;
            else
                frc.Shortened = 0;
            end

            % Set up variable indicating the current transport block sizes
            TBSs = trBlkSizes(:, mod(subframeNo, 10)+1).';

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

            % Equalization and combining, layer demapping, transform
            % deprecoding, demodulation and descrambling of the received
            % data
            ulschDims = lteULSCHInfo(frc, frc.PUSCH, TBSs, 'chsconcat');
            if (strcmpi(reference, 'Antennas')==1)
                % MMSE equalization occurs inside ltePUSCHDecode
                rxEncodedBits = ltePUSCHDecode(frc, ulschDims, ...
                    puschRx, puschEstCh, noiseEst);
            elseif (strcmpi(reference, 'Layers')==1)
                rxSymbols = lteEqualizeMMSE(puschRx, puschEstCh, noiseEst);
                rxEncodedBits = ltePUSCHDecode(frc, ulschDims, rxSymbols);
            end

            % UL-SCH transport decoding
            [rxdata, harqProcesses(harqID).blkerr, ...
                harqProcesses(harqID).decState] = lteULSCHDecode( ...
                frc, ulschDims, TBSs, rxEncodedBits, ...
                harqProcesses(harqID).decState);

            % Store block CRC and throughput results for subframes
            % containing transport data
            crcSNR(dataSubframeIndex, :) = harqProcesses(harqID).blkerr;
            tputSNR(dataSubframeIndex, :) = TBSs.* ...
                (1-harqProcesses(harqID).blkerr);
            dataSubframeIndex = dataSubframeIndex + 1;

            % Provide PMI feedback (if multiple antennas are being used)
            if (frc.NTxAnts>1)
                if (~any(TBSs) || ~updatePMI)
                    % Use previous PMI value if in a subframe not used for
                    % uplink in TDD, or if PMI update is not configured for
                    % this subframe.
                    PMI = txPMIs(:, mod(pmiIdx-1, pmiDelay)+1);
                else
                    % Update PMI estimate. If not using the perfect channel
                    % estimator, redo channel estimation using only the SRS
                    % for PMI selection purposes.
                    if (~perfectChanEstimator)
                        cec.Reference = 'None';
                        [estChannelGrid, noiseEst] = ...
                            lteULChannelEstimate(frc, frc.PUSCH, cec, ...
                            rxSubframe, refGrid);
                    end
                    PMI = lteULPMISelect(frc, frc.PUSCH, ...
                        estChannelGrid, noiseEst, refGrid, cec);
                end
                txPMIs(:, pmiIdx+1) = PMI;
            end
        end
    end

    % Record the block CRC error and bit throughput for the total number of
    % frames simulated at an SNR point
    crc(ResultIndex, :) = crcSNR(:);
    tput(ResultIndex, :) = tputSNR(:);
    ResultIndex = ResultIndex + 1;
    disp(' ');

end

hMultiCodewordPUSCHResults(SNRIn, NFrames, trBlkSizes, crc, tput);