NFrames = 10;
SNRdB = 6.0;

enb = struct('RC','R.3');         % Set up parameters of RMC R.3
enb = lteRMCDL(enb);
enb.CFI = 3;                      % Reconfigure Control Format Indicator
enb.OCNGPDSCHEnable = 'On';       % Enable OCNG for unallocated PDSCH REs
enb.TotSubframes = 1;             % Reconfigure for a single subframe
enb.PDSCH.RVSeq = 0;              % Disable HARQ
enb.PDSCH.CSIMode = 'PUCCH 1-0';  % Configure the CSI reporting mode
enb.PDSCH.CSI = 'On';             % CSI scaling of soft bits

channel.Seed = 10;                  % Channel seed
channel.NRxAnts = 2;                % 2 receive antennas
channel.DelayProfile = 'EPA';       % Delay profile
channel.DopplerFreq = 5.0;          % Doppler frequency
channel.MIMOCorrelation = 'High';   % Multi-antenna correlation
channel.ModelType = 'GMEDS';        % Rayleigh fading model type
channel.NormalizeTxAnts = 'On';     % Normalize for transmit antennas
channel.NormalizePathGains = 'On';  % Normalize delay profile power
channel.InitPhase = 'Random';       % Random initial phases
channel.NTerms = 16;                % Oscillators used in fading model

 % Set channel model sampling rate
ofdmInfo = lteOFDMInfo(enb);
channel.SamplingRate = ofdmInfo.SamplingRate;

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

% Channel estimator behavior
perfectChanEstimator = true;

cqiDelay = 8; % subframes

% Initialize variables used for results recording
CQIReport = [];     % reported CQI values
SINRReport = [];    % corresponding SINR values
xaxis = [];         % corresponding subframe numbers

% For each CQI configuration (UE reported and median):
for cqiConfig = 1:2

    if (cqiConfig==1)
        cqiConfigStr = 'UE reported';
    else
        cqiConfigStr = 'median';
    end
    fprintf('\nSimulating with %s CQI at %gdB SNR for %d Frame(s)\n', ...
        cqiConfigStr,SNRdB,NFrames);

    % Initialize CQI values: for UE reported, set to all ones; for median,
    % set to the median of the CQI values for the UE reported run
    if (cqiConfig==1)
        cqiBuffer = ones(1,cqiDelay);
    else
        cqiBuffer = ones(1,cqiDelay)*medianCQI;
    end

    % Initialize variables
    rng('default');     % Default random number generator seed
    totalCRC = [];      % CRC values, used for throughput calculation
    totalTBS = [];      % TBS values, used for throughput calculation
    offsets = 0;        % Initialize frame offset value

    % For each subframe:
    for subframeNo = 0:(NFrames*10-1)

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

        % Select CQI, reading the oldest value from the CQI buffer
        cqiPtr = mod(subframeNo,cqiDelay);
        CQI = cqiBuffer(cqiPtr+1);

        % Select MCS according to CQI using TS36.101 Table A.4-1 CSI RMC
        % RC.1 FDD (MCS.1), which defines the relationship between CQI
        % indices and MCS indices
        IMCSTable = [-1 0 0 2 4 6 8 11 13 16 18 21 23 25 27 27];
        IMCS = IMCSTable(CQI+1);

        % Determine TBS and modulation order
        [ITBS,modulation] = lteMCS(IMCS);
        enb.PDSCH.Modulation = {modulation};
        if (mod(enb.NSubframe,5)==0)
            TBS = 0;
        else
            TBS = double(lteTBS(size(enb.PDSCH.PRBSet,1),ITBS));
        end
        enb.PDSCH.TrBlkSizes(enb.NSubframe+1) = TBS;

        % Determine if a CQI update is required in this subframe, according
        % to reporting periodicity N_pd = 2ms and configuration index
        % cqi-pmi-ConfigurationIndex = 1 from TS36.101 Table 9.3.2.1.1-1
        cqiPeriod = 2; % periodicity N_pd
        cqiOffset = 1; % offset deriving from cqi-pmi-ConfigurationIndex
        cqiUpdate = (mod(subframeNo,cqiPeriod)==cqiOffset);

        % Establish if this subframe actually needs executed for PDSCH
        % reception, CQI estimation or initial timing offset estimation:
        if((TBS~=0 && subframeNo>=(cqiDelay+cqiOffset)) || ...
                (cqiConfig==1 && cqiUpdate) || subframeNo==0)

            % Generate random bits for the subframe
            data = randi([0 1],TBS,1);

            % Create OFDM resource grid containing RMC transmission and
            % perform OFDM modulation.
            txWaveform = lteRMCDLTool(enb,data);

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

            % Pass data through the fading channel model.
            % An additional 25 samples are added to the end of the
            % waveform. These are to cover the range of delays expected
            % from the channel modeling (a combination of implementation
            % delay and channel delay spread)
            rxWaveform = lteFadingChannel(channel, ...
                            [txWaveform ; zeros(25,size(txWaveform,2))]);

            % 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(ofdmInfo.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;

            % Perform synchronization
            % 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)
                offset = lteDLFrameOffset(enb,rxWaveform);
                if (offset > 25)
                    offset = offsets(end);
                end
                offsets = [offsets offset]; %#ok<AGROW>
            end
            rxWaveform = rxWaveform(1+offset:end,:);

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

            % Channel estimation
            if (perfectChanEstimator)
                chEstGrid = ...
                   lteDLPerfectChannelEstimate(enb,channel,offset);
                n = lteOFDMDemodulate(enb,noise(1+offset:end,:));
                noiseEst = var(reshape(n,numel(n),1));
            else
                [chEstGrid,noiseEst] = ...
                    lteDLChannelEstimate(enb,enb.PDSCH, ...
                        cec,rxGrid); %#ok<UNRCH>
            end

            % If this subframe requires PDSCH reception:
            if(TBS~=0 && subframeNo>=(cqiDelay+cqiOffset))

                % Decode the PDSCH
                ind = ltePDSCHIndices(enb,enb.PDSCH,enb.PDSCH.PRBSet);
                pdschRx = lteExtractResources(ind,rxGrid) * ...
                        (10^(-enb.PDSCH.Rho/20));
                pdschChEst = lteExtractResources(ind,chEstGrid);
                [rxBits,rxSymbols] = ltePDSCHDecode(enb,enb.PDSCH, ...
                                        pdschRx,pdschChEst,noiseEst);

                % Decode the DL-SCH
                [decbits,crc] = lteDLSCHDecode(enb,enb.PDSCH,TBS,rxBits);

                % Record the CRC and TBS values for final throughput
                % calculation
                totalCRC = [totalCRC crc]; %#ok<AGROW>
                totalTBS = [totalTBS TBS]; %#ok<AGROW>

            end

            % Update CQI:
            if (cqiConfig==1 && cqiUpdate)

                % Perform CQI selection
                [thisCQI,thisSINR] = ...
                    lteCQISelect(enb,enb.PDSCH,chEstGrid,noiseEst);

                % Feed the CQI value back to UE (in a buffer)
                cqiBuffer(cqiPtr+1) = thisCQI;

                % Record values for plotting
                CQIReport = [CQIReport thisCQI];    %#ok<AGROW>
                SINRReport = [SINRReport thisSINR]; %#ok<AGROW>
                xaxis = [xaxis subframeNo];         %#ok<AGROW>

            end

        end

        % For subframes where CQI was not updated, re-use the previous
        % value in the buffer
        if (cqiConfig==1 && ~cqiUpdate)
            cqiBuffer(cqiPtr+1) = cqiBuffer(mod(cqiPtr-1,cqiDelay)+1);
        end

    end

    % Display results for the current CQI configuration
    fprintf('\nResults with %s CQI:\n',cqiConfigStr);
    tputTotal = sum(totalTBS);
    if (cqiConfig==1)

        % Compute and display throughput
        tputUEReported = sum(totalTBS.*(1-totalCRC));
        fprintf('Throughput: %d bits (%0.2f%%)\n', ...
            tputUEReported,tputUEReported/tputTotal*100);

        % Compute and display BLER
        measuredBLER = mean(totalCRC);
        fprintf('BLER: %0.3f (requirement is >= 0.02)\n',measuredBLER);

        % Compute and display median CQI
        medianCQI = ceil(median(CQIReport));
        fprintf('Median CQI: %d\n',medianCQI');

        % Compute and display proportion of CQI values
        % outside +/- 1 of the median
        measuredAlpha = (sum(CQIReport<(medianCQI-1)) + ...
             sum(CQIReport>(medianCQI+1)))/length(CQIReport)*100;
        fprintf(['Percentage of CQI indices outside +/- 1 of median:' ...
            ' %0.2f%% (requirement is >= 20%%)\n'],measuredAlpha);

    else

        % Compute and display throughput
        tputMedian = sum(totalTBS.*(1-totalCRC));
        fprintf('Throughput: %d bits (%0.2f%%)\n', ...
            tputMedian,tputMedian/tputTotal*100);

        % Compute and display throughput ratio
        measuredGamma = tputUEReported/tputMedian;
        fprintf(['Throughput ratio (gamma): %0.3f' ...
            ' (requirement is >= 1.05)'],measuredGamma);

    end

end

figure;
subplot(2,1,1);
plot(xaxis,SINRReport,'rx-');
xlabel('Subframe number');
ylabel('Estimated SINR (dB)');
title('Estimated SINR versus subframe number');
hold on;
subplot(2,1,2);
plot(xaxis,CQIReport,'bo-');
xlabel('Subframe number');
ylabel('Selected CQI Index');
title('Selected CQI Index versus subframe number');