NFrames = 10;
SNRdB = 0.0;
FixedRI = 2;

channel.Seed = 1;                   % Channel seed
channel.NRxAnts = 2;                % 2 receive antennas
channel.DelayProfile = 'EPA';       % Delay profile
channel.DopplerFreq = 5.0;          % Doppler frequency
channel.MIMOCorrelation = 'Low';    % 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

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

riDelay = 9; % subframes
cqipmiDelay = 8; % subframes

codebookSubsetRI1 = '000011';
codebookSubsetRI2 = '010000';
codebookSubsetUEReported = '010011';

% Check that fixed RI valid is either 1 or 2
if (~any(FixedRI==[1 2]))
    error('Fixed RI value specified must be 1 or 2.');
end

% For each RI configuration (fixed RI and UE reported):
for riConfig = 1:2

    if (riConfig==1)
        riConfigStr = sprintf('fixed RI=%d',FixedRI);
    else
        riConfigStr = 'UE reported RI';
    end
    fprintf('\nSimulating with %s at %gdB SNR for %d Frame(s)\n', ...
        riConfigStr,SNRdB,NFrames);

    % Set up eNodeB settings for 2 codewords
    enb = struct('RC','R.3');
    enb.CellRefP = 2;
    enb.CFI = 3;
    enb.OCNGPDSCHEnable = 'On';
    enb.TotSubframes = 1;
    enb.PDSCH.RVSeq = 0;
    enb.PDSCH.CSIMode = 'PUCCH 1-1';
    enb.PDSCH.TxScheme = 'SpatialMux';
    enb.PDSCH.NLayers = 2;
    enb.PDSCH.Rho = -3.0;
    enb.PDSCH.CSI = 'On';
    ncw = 2;
    enb = lteRMCDL(enb,ncw);

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

    % Configure appropriate codebook subset restriction and target
    % throughput ratio
    if (riConfig==1)
        if (FixedRI==1)
            enb.PDSCH.CodebookSubset = codebookSubsetRI1;
            targetGamma = 1.05;
        else
            enb.PDSCH.CodebookSubset = codebookSubsetRI2;
            targetGamma = 1.0;
        end
    else
        enb.PDSCH.CodebookSubset = codebookSubsetUEReported;
    end

    % Initialize variables
    totalCRC = [];      % CRC values, used for throughput calculation
    totalTBS = [];      % TBS values, used for throughput calculation
    RIReport = [];      % reported RI values
    riXaxis = [];       % corresponding subframe numbers
    PMIReport = [];     % reported PMI values
    CQIReport = [];     % reported CQI values
    SINRReport = [];    % corresponding SINR values
    cqipmiXaxis = [];   % corresponding subframe numbers
    offsets = 0;        % Initialize frame offset value
    rng('default');     % Default random number generator seed

    % Initialize RI/PMI/CQI buffers; note that the PDSCH throughput
    % resulting from these initial values is ignored, as the throughput
    % recording waits until the buffers have been filled with meaningful
    % reports
    riBuffer = ones(1,riDelay);
    pmiBuffer = ones(1,cqipmiDelay);
    cqiBuffer = ones(ncw,cqipmiDelay);

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

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

        % Update PMI, reading the oldest value from the PMI buffer
        cqipmiPtr = mod(subframeNo,cqipmiDelay);
        enb.PDSCH.PMISet = pmiBuffer(cqipmiPtr+1);

        % Update number of transmission layers based on RI, reading the
        % oldest value from the RI buffer
        riPtr = mod(subframeNo,riDelay);
        enb.PDSCH.NLayers = riBuffer(riPtr+1);

        % Update MCS according to CQI, reading the oldest value from the
        % CQI buffer
        CQI = cqiBuffer(:,cqipmiPtr+1);
        % Select MCS according to CQI using TS36.101 Table A.4-1 CSI RMC
        % RC.2 FDD (MCS.2), which defines the relationship between CQI
        % indices and MCS indices
        IMCSTable = [-1 0 0 2 4 6 8 11 13 15 18 20 22 24 26 27];
        IMCS = IMCSTable(CQI+1);

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

        % Determine if an RI or PMI/CQI update is required in this
        % subframe, according to reporting periodicity N_pd = 5ms and
        % configuration indices cqi-pmi-ConfigurationIndex = 6 and
        % ri-ConfigurationInd = 1 from TS36.101 Table 9.5.1.1-1
        riPeriod = 5; % periodicity N_pd
        riOffset = 3; % offset deriving from configuration indices
        riUpdate = (mod(subframeNo,riPeriod)==riOffset);
        cqipmiPeriod = 5; % periodicity N_pd
        cqipmiOffset = 4; % offset deriving from cqi-pmi-ConfigurationIndex
        cqipmiUpdate = (mod(subframeNo,cqipmiPeriod)==cqipmiOffset);

        % Establish if this subframe actually needs executed for PDSCH
        % reception, RI/PMI/CQI estimation or initial timing offset
        % estimation:
        firstRxSubframe = max(cqipmiDelay+cqipmiOffset,riDelay+riOffset);
        if((any(TBS) && subframeNo>=firstRxSubframe) || ...
                (riUpdate || cqipmiUpdate) || subframeNo==0)

            % Configure data for transmission
            % (1 codeword on 1 layer, 2 codewords on 2 layers)
            data = {randi([0 1],TBS(1),1)};
            if (enb.PDSCH.NLayers==2)
                data = {data{1} randi([0 1],TBS(2),1)};
            end
            enb.PDSCH.Modulation = modulation(1:numel(data));

            % 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(any(TBS) && subframeNo>=firstRxSubframe)

                % 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, ...
                    enb.PDSCH.TrBlkSizes(1:numel(rxBits),enb.NSubframe+1),rxBits);

                % Record the CRC and TBS values for final throughput
                % calculation
                totalCRC = [totalCRC crc];                  %#ok<AGROW>
                totalTBS = [totalTBS TBS(1:numel(crc)).'];  %#ok<AGROW>

            end

            % Update RI:
            if (riUpdate)

                % Update RI
                thisRI = lteRISelect(enb,enb.PDSCH,chEstGrid,noiseEst);

                % Feed the value back to UE (in a buffer)
                riBuffer(riPtr+1) = thisRI;

                % Record values for plotting
                if (riConfig==2)
                    RIReport = [RIReport thisRI];           %#ok<AGROW>
                    riXaxis = [riXaxis subframeNo];         %#ok<AGROW>
                end

            end

            % Update CQI and PMI:
            if (cqipmiUpdate)

                % Update PMI conditioned on the most recently reported RI,
                % or the fixed RI if RI reporting is not configured
                if (isempty(RIReport))
                    enb.PDSCH.NLayers = FixedRI;
                else
                    enb.PDSCH.NLayers = RIReport(end);
                end
                thisPMI = ltePMISelect(enb,enb.PDSCH,chEstGrid,noiseEst);

                % Update CQI conditioned on the most recently reported PMI
                % and RI; this includes configuring the number of codewords
                % for CQI selection based on the RI
                enb.PDSCH.PMISet = thisPMI;
                enb.PDSCH.NCodewords = enb.PDSCH.NLayers;
                [thisCQI,thisSINR] = ...
                    lteCQISelect(enb,enb.PDSCH,chEstGrid,noiseEst);

                % Feed the values back to UE (in buffers)
                pmiBuffer(:,cqipmiPtr+1) = thisPMI;
                % For CSIMode='PUCCH 1-1', the CQI for the second codeword
                % is reported as a differential from the first codeword;
                % here we convert to absolute CQI indices for both
                % codewords.
                if (numel(thisCQI)==2)
                    thisCQI(2) = thisCQI(1) - thisCQI(2);
                end
                % Transmit with lowest CQI=1 when CQI selection reports
                % CQI=0 (out of range)
                thisCQI(thisCQI==0) = 1;

                cqiBuffer(:,cqipmiPtr+1) = thisCQI;

                % Record values for plotting
                if (riConfig==2)
                    PMIReport = [PMIReport thisPMI];        %#ok<AGROW>
                    CQIReport = [CQIReport thisCQI(1)];     %#ok<AGROW>
                    SINRReport = [SINRReport thisSINR(1)];  %#ok<AGROW>
                    cqipmiXaxis = [cqipmiXaxis subframeNo]; %#ok<AGROW>
                end

            end

        end

        % For subframes where RI was not updated, re-use the previous value
        % in the buffer
        if (~riUpdate)
            riBuffer(riPtr+1) = riBuffer(mod(riPtr-1,riDelay)+1);
        end

        % For subframes where PMI/CQI were not updated, re-use the previous
        % values in the buffers
        if (~cqipmiUpdate)
            pmiBuffer(:,cqipmiPtr+1) = ...
                pmiBuffer(:,mod(cqipmiPtr-1,cqipmiDelay)+1);
            cqiBuffer(:,cqipmiPtr+1) = ...
                cqiBuffer(:,mod(cqipmiPtr-1,cqipmiDelay)+1);
        end

    end

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

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

    else

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

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

    end

end

figure;
subplot(2,1,1);
plot(riXaxis,RIReport,'bo-');
xlabel('Subframe number');
ylabel('Selected RI');
title('Selected RI versus subframe number');
hold on;
subplot(2,1,2);
plot(cqipmiXaxis,PMIReport,'ks-');
xlabel('Subframe number');
ylabel('Selected PMI');
title('Selected PMI versus subframe number');

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