NFrames = 4;               % Number of frames
eqMethod = 'MMSE_IRC';     % MMSE, MMSE_IRC

SINR = 0.8;        % SINR in dB
DIP2 = -1.73;      % DIP in dB for cell 2
DIP3 = -8.66;      % DIP in dB for cell 3
Noc = -98;         % dBm/15kHz average power spectral density

% Set the random number generator seed
rng('default');

% Set cell 1 eNodeB configuration according to R.47
simulationParameters = struct;
simulationParameters.NDLRB = 50;
simulationParameters.CellRefP = 2;
simulationParameters.NCellID = 0;
simulationParameters.CFI = 2;
simulationParameters.DuplexMode = 'FDD';
simulationParameters.TotSubframes = 1; % This is not the total number of
% subframes used in the simulation, just the number of subframes we
% generate per call to the waveform generator.

% Specify PDSCH configuration substructure
simulationParameters.PDSCH.TxScheme = 'SpatialMux';
simulationParameters.PDSCH.Modulation = {'16QAM'};
simulationParameters.PDSCH.NLayers = 1;
simulationParameters.PDSCH.Rho = -3;
simulationParameters.PDSCH.PRBSet = (0:49)';
% Table A.3.3.2.1-2, TS 36.101
simulationParameters.PDSCH.TrBlkSizes = [8760 8760 8760 8760 8760 0 8760 8760 8760 8760];
% Table A.3.3.2.1-2, TS 36.101
simulationParameters.PDSCH.CodedTrBlkSizes = [24768 26400 26400 26400 26400 0 26400 26400 26400 26400];
simulationParameters.PDSCH.CSIMode = 'PUCCH 1-1';
simulationParameters.PDSCH.PMIMode = 'Wideband';
simulationParameters.PDSCH.CSI = 'On';
simulationParameters.PDSCH.W = [];
simulationParameters.PDSCH.CodebookSubset = '1111';

% Specify PDSCH OCNG configuration
simulationParameters.OCNGPDSCHEnable = 'On';             % Enable OCNG fill
simulationParameters.OCNGPDSCHPower = -3;                % OCNG power same as PDSCH Rho
simulationParameters.OCNGPDSCH.RNTI = 0;                 % Virtual UE RNTI
simulationParameters.OCNGPDSCH.Modulation = 'QPSK';      % OCNG symbol modulation
simulationParameters.OCNGPDSCH.TxScheme = 'TxDiversity'; % OCNG transmission mode 2

enb1 = lteRMCDL(simulationParameters);

% Cell 2
enb2 = enb1;
enb2.NCellID = 1;
enb2.OCNGPDSCHEnable = 'Off';

% Cell 3
enb3 = enb1;
enb3.NCellID = 2;
enb3.OCNGPDSCHEnable = 'Off';

% eNodeB1 to UE propagation channel
channel1 = struct;                    % Channel config structure
channel1.Seed = 20;                   % Channel seed
channel1.NRxAnts = 2;                 % 2 receive antennas
channel1.DelayProfile ='EVA';         % Delay profile
channel1.DopplerFreq = 5;             % Doppler frequency
channel1.MIMOCorrelation = 'Low';     % Multi-antenna correlation
channel1.NTerms = 16;                 % Oscillators used in fading model
channel1.ModelType = 'GMEDS';         % Rayleigh fading model type
channel1.InitPhase = 'Random';        % Random initial phases
channel1.NormalizePathGains = 'On';   % Normalize delay profile power
channel1.NormalizeTxAnts = 'On';      % Normalize for transmit antennas

ofdmInfo = lteOFDMInfo(enb1);
channel1.SamplingRate = ofdmInfo.SamplingRate;

% eNodeB2 (interference) to UE propagation channel
channel2 = channel1;
channel2.Seed = 122;                  % Channel seed

% eNodeB3 (interference) to UE propagation channel
channel3 = channel1;
channel3.Seed = 36;                   % Channel seed

% Channel estimator behavior
perfectChanEstimator = true;

% The channel estimator configuration structure is defined below. The
% frequency and time averaging window sizes are chosen to span a relatively
% large number of resource elements. The large window sizes are chosen to
% average as much as possible noise and interference in the resource
% elements. Note that too large an averaging window in time and/or frequency
% will cause loss of information due to averaging out the channel variations.
% This produces an increasingly imperfect channel estimate which can affect
% the performance of the equalizer.

cec = struct;                        % Channel estimation config structure
cec.PilotAverage = 'UserDefined';    % Type of pilot symbol averaging
cec.FreqWindow = 31;                 % Frequency window size
cec.TimeWindow = 23;                 % Time window size
cec.InterpType = 'Cubic';            % 2D interpolation type
cec.InterpWindow = 'Centered';       % Interpolation window type
cec.InterpWinSize = 1;               % Interpolation window size

% Channel noise setup
nocLin = 10.^(Noc/10)*(1e-3); % linear  in Watts
% Take into account FFT (OFDM) scaling
No = sqrt(nocLin/(2*double(ofdmInfo.Nfft)));

% Signal and interference amplitude scaling factors calculation
[K1, K2, K3] = hENBscalingFactors(DIP2, DIP3, Noc, SINR, enb1, enb2, enb3);

% Initialize state of all HARQ processes
harqProcesses = hNewHARQProcess(enb1);
% 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;

% Set up variables for the main loop
lastOffset = 0;       % Initialize frame timing offset from previous frame
frameOffset = 0;      % Initialize frame timing offset
nPMI = 0;             % Initialize the number of precoder matrix indication
                      % (PMI) set calculated
blkCRC = [];          % Block CRC for all considered subframes
bitTput = [];         % Number of successfully received bits per subframe
txedTrBlkSizes = [];  % Number of transmitted bits per subframes
% Vector of total number of bits transmitted, calculated for each subframe.
runningMaxThPut = [];
% Vector storing the number of successfully received bits, calculated for
% each subframe.
runningSimThPut = [];

% Obtain the number of transmit antennas.
dims = lteDLResourceGridSize(enb1);
P = dims(3);

% Assign the redundancy version sequence for each codeword and transport
% block sizes for each subframe
rvSequence = enb1.PDSCH.RVSeq;
trBlkSizes = enb1.PDSCH.TrBlkSizes;

% Set the PMI delay for the closed-loop spatial multiplexing
pmiDelay = 8; % As specified in TS 36.101, Table 8.2.1.4.1B-1

% Initialize PMIs for the first 'pmiDelay' subframes
pmiDims = ltePMIInfo(enb1,enb1.PDSCH);
txPMIs = zeros(pmiDims.NSubbands, pmiDelay);

% Flag to indicate if a valid PMI feedback is available from the UE
pmiReady = false;

fprintf('\nSimulating %d frame(s)\n',NFrames);

% Main for loop: for all subframes
for subframeNo = 0:(NFrames*10-1)

    % Reinitialize channel seed for each subframe to increase variability
    channel1.Seed = 1+subframeNo;
    channel2.Seed = 1+subframeNo+(NFrames*10);
    channel3.Seed = 1+subframeNo+2*(NFrames*10);

    % Update subframe number
    enb1.NSubframe = subframeNo;
    enb2.NSubframe = subframeNo;
    enb3.NSubframe = subframeNo;

    duplexInfo = lteDuplexingInfo(enb1);

    if  duplexInfo.NSymbolsDL ~= 0 % target only downlink subframes

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

        % Update the PDSCH transmission config with HARQ process state
        enb1.PDSCH.RVSeq = harqProcesses(harqID).txConfig.RVSeq;
        enb1.PDSCH.RV = harqProcesses(harqID).txConfig.RV;
        dlschTransportBlk = harqProcesses(harqID).data;

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

        % Create transmit waveform
        [tx,~,enbOut] = lteRMCDLTool(enb1, dlschTransportBlk);

        % Pad 25 samples to cover the range of delays expected from channel
        % modeling (a combination of implementation delay and channel delay
        % spread)
        txWaveform1 = [tx; zeros(25, P)];

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

        % Generate interferer model as per as per TS 36.101, B.5.3. The
        % function hTM4InterfModel generates the interferer transmit signal.
        txWaveform2 = [hTM4InterfModel(enb2); zeros(25,P)];
        txWaveform3 = [hTM4InterfModel(enb3); zeros(25,P)];

        % Specify channel time for the present subframe
        channel1.InitTime = subframeNo/1000;
        channel2.InitTime = channel1.InitTime;
        channel3.InitTime = channel1.InitTime;

        % Pass data through the channel
        rxWaveform1 = lteFadingChannel(channel1,txWaveform1);
        rxWaveform2 = lteFadingChannel(channel2,txWaveform2);
        rxWaveform3 = lteFadingChannel(channel3,txWaveform3);

        % Generate noise
        noise = No*complex(randn(size(rxWaveform1)), ...
            randn(size(rxWaveform1)));

        % Add AWGN to the received time domain waveform
        rxWaveform = K1*rxWaveform1 + K2*rxWaveform2 + K3*rxWaveform3 + noise;

        % Receiver
        % Once every frame, on subframe 0, calculate a new synchronization
        % offset
        if (mod(subframeNo,10) == 0)
            frameOffset = lteDLFrameOffset(enb1, rxWaveform);
            if (frameOffset > 25)
                frameOffset = lastOffset;
            end
            lastOffset = frameOffset;
        end

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

        % Scale rxWaveform by 1/K1 to avoid numerical issues with
        % channel decoding stages
        rxWaveform = rxWaveform/K1;

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

        % Perform channel estimation
        if(perfectChanEstimator)
            estChannelGrid = lteDLPerfectChannelEstimate(enb1, channel1, frameOffset);
            noiseInterf = K2*rxWaveform2 + K3*rxWaveform3 + noise;
            noiseInterf = noiseInterf/K1;
            noiseGrid = lteOFDMDemodulate(enb1, noiseInterf(1+frameOffset:end ,:));
            noiseEst = var(noiseGrid(:));
        else
            [estChannelGrid, noiseEst] = lteDLChannelEstimate( ...
                enb1, enb1.PDSCH, cec, rxSubframe);
        end

        % Get PDSCH indices
        pdschIndices = ltePDSCHIndices(enb1,enb1.PDSCH,enb1.PDSCH.PRBSet);
        % Get PDSCH resource elements. Scale the received subframe by
        % the PDSCH power factor Rho.
        [pdschRx, pdschHest] = lteExtractResources(pdschIndices, ...
            rxSubframe*(10^(-enb1.PDSCH.Rho/20)), estChannelGrid);

        % Perform equalization and deprecoding
        if strcmp(eqMethod,'MMSE')
            % MIMO equalization and deprecoding (MMSE based)
            [rxDeprecoded,csi] = lteEqualizeMIMO(enb1,enb1.PDSCH,...
                pdschRx,pdschHest,noiseEst);
        else
            % MIMO equalization and deprecoding (MMSE-IRC based)
            [rxDeprecoded,csi] = hEqualizeMMSEIRC(enb1,enb1.PDSCH,...
                rxSubframe,estChannelGrid,noiseEst);
        end

        % Perform layer demapping, demodulation and descrambling
        cws = ltePDSCHDecode(enb1,setfield(enb1.PDSCH,'TxScheme',...
            'Port7-14'),rxDeprecoded); % The PDSCH transmission scheme is
        % modified into port7-14, in order to skip the deprecoding operation

        % Scaling LLRs by CSI
        cws = hCSIscaling(enb1.PDSCH,cws,csi);

        % Decode DL-SCH
        [decbits, harqProcesses(harqID).blkerr,harqProcesses(harqID).decState] = ...
            lteDLSCHDecode(enb1, enb1.PDSCH, trBlk, cws, ...
            harqProcesses(harqID).decState);

        % Store values to calculate throughput
        % Only for subframes with data and valid PMI feedback
        if any(trBlk) && pmiReady
            blkCRC = [blkCRC harqProcesses(harqID).blkerr];
            txedTrBlkSizes =  [txedTrBlkSizes trBlk];
            bitTput = [bitTput trBlk.*(1-harqProcesses(harqID).blkerr)];
        end
        runningSimThPut = [runningSimThPut sum(bitTput,2)];
        runningMaxThPut = [runningMaxThPut sum(txedTrBlkSizes,2)];

        % Provide PMI feedback to the eNodeB
        if strcmpi(enb1.PDSCH.TxScheme,'SpatialMux')
            PMI = ltePMISelect(enb1, enb1.PDSCH, estChannelGrid, noiseEst);
            txPMIs(:, pmiIdx+1) = PMI;
            nPMI = nPMI+1;
            if nPMI>=pmiDelay
                pmiReady = true;
            end
        end
    end
end

maxThroughput = sum(txedTrBlkSizes); % Maximum possible throughput
simThroughput = sum(bitTput,2);      % Simulated throughput

% Display achieved throughput percentage
disp(['Achieved throughput ' num2str(simThroughput*100/maxThroughput) '%'])

% Plot running throughput
figure;plot(runningSimThPut*100./runningMaxThPut)
ylabel('Throughput (%)');
xlabel('Simulated subframe');
title('Throughput');