% eNodeB Configuration
enb = struct;                       % eNodeB config structure
enb.RC = 'R.12';                    % RMC number
enb.NCellID = 10;                   % Cell ID
% Total number of subframes lteRMCDLTool  will generate per call. This
% example generates 1 frame (10 subframes) worth of data each call within
% the processing loop. lteRMCDLTool will be called 10 (NFrames) times.
enb.TotSubframes = 10;              
enb.PDSCH.TxScheme = 'TxDiversity'; % Transmission scheme
enb.PDSCH.RNTI = 1;                 % 16-bit UE specific mask
enb.PDSCH.Rho = -3;                 % Downlink power allocation 
enb.PDSCH.CSI = 'On';               % CSI scaling of soft bits
enb.PDSCH.RVSeq = 0;                % Disable HARQ

NFrames = 10;                       % Number of frames to simulate
SNRIn = [-2 2];                     % SNR points to simulate

% Generate RMC configuration for test model specified in enb structure
rmc = lteRMCDL(enb);
    
% Transport block sizes of each subframe within a frame
trblksize = rmc.PDSCH.TrBlkSizes;
ncw = size(trblksize,1);             % no of codewords associated to RMC

% dims is a three element vector [K; L; P]: where K = no of subcarriers,
% L = no of OFDM symbols and P = no of transmit antennas
dims = lteDLResourceGridSize(rmc);  % Determine resource grid dimensions
L = dims(2);                        % Number of OFDM symbols per subframe
P = dims(3);                        % Number of transmit antennas

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

cec = struct;                       % Channel estimation config structure
cec.PilotAverage = 'UserDefined';   % Type of pilot symbol averaging
cec.FreqWindow = 15;                % Frequency window size in REs
cec.TimeWindow = 141;               % Time window size in REs

cec.InterpType = 'Cubic';         % Cubic interpolation
cec.InterpWinSize = 3;            % Interpolate up to 3 subframes 
                                  % simultaneously
cec.InterpWindow = 'Centred';     % Interpolation windowing method

% Initialize the variables used in the simulation and analysis
resultIndex = 1;                    % Initialize SNR counter index
offsets = 0;                        % Initialize offset vector
crc = [];                           % Define CRC error per frame
trSizes = [];                       % Define transport sizes per frame
    
% Set the random number generator to default value
rng('default');

% Number of subframes carrying transport blocks with data within a frame.
% For RMC R.12 this should be 9, of the 10 subframes in a frame subframe 5
% does not carry any data
nDataTBSPerFrame = sum(trblksize(:) ~= 0);

% Total block CRC vector
TotalBLKCRC = zeros(numel(SNRIn),nDataTBSPerFrame*NFrames); 

% Bit throughput vector
BitThroughput = zeros(numel(SNRIn),nDataTBSPerFrame*NFrames);

for SNRdB = SNRIn
    fprintf('\nSimulating at %gdB SNR for a total %d Frame(s)\n',...
    SNRdB,NFrames);

    % Initialize result store for each SNR point
    TotalBLKCRCcws = zeros(NFrames,nDataTBSPerFrame);% Total block CRC 
                                                     % vector
    BitTput = zeros(NFrames,nDataTBSPerFrame);      % Bit throughput vector

    for FrameNo = 1:NFrames
        % Generate random bits for the frame. The number of bits to
        % generate for each subframe is specified by each element of
        % trblksize. For a full frame generate sum(trblksize) bits
        trdata = randi([0 1], sum(trblksize), 1);
        
        % Generate populated LTE resource grid using RMC generator and OFDM
        % modulate. info is a structure containing the sampling rate used
        % for OFDM modulation
        [txWaveform,txGrid,info] = lteRMCDLTool(rmc,trdata);
       
        % Initialize result store for frame
        crc = zeros(nDataTBSPerFrame/ncw,ncw);    % Intermediate block CRC 
        trSizes = zeros(nDataTBSPerFrame/ncw,ncw);% Intermediate throughput
        dataSubframeIndex = 1;
        
        % Set sampling rate of channel to that of OFDM modulation     
        cfg.SamplingRate = info.SamplingRate;

        % Set channel offset to current frame (1 frame = 10ms)
        cfg.InitTime = (FrameNo-1)*(rmc.TotSubframes)/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(cfg,[txWaveform ; zeros(25,P)]);

        % Noise setup
        SNR = 10^(SNRdB/20);    % Linear SNR
        
        % 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*rmc.CellRefP*double(info.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 receiver synchronization 
        offset = lteDLFrameOffset(rmc,rxWaveform);

        % Determine if frame synchronization was successful
        if (offset > 25)
            offset = offsets(end);
        else
            offsets = [offsets offset]; %#ok
        end
        if (offset>0)
            rxWaveform = rxWaveform(1+offset:end,:);
        end

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

        % Perform channel estimation
        [estChannelGrid, noiseest] = lteDLChannelEstimate(rmc,cec,rxGrid);
        
        % Process subframes 0 to 9 within received frame
        for sf = 0:rmc.TotSubframes-1
            
            % Increment NSubframe for correct decoding of data
            rmc.NSubframe = mod(sf,10);

            % Extract one subframe for analyzing at a time from the
            % received grid        
            rxSubframe = rxGrid(:,L*sf+1:L*(sf+1),:);
            
            % Extract the estimated channel response for the subframe
            % being analyzed
            chSubframe = estChannelGrid(:,L*sf+1:L*(sf+1),:,:);

            % Perform deprecoding, layer demapping, demodulation and
            % descrambling on the received data using the estimate of
            % the channel
            rxEncodedBits = ltePDSCHDecode(rmc,rmc.PDSCH,...
                  rxSubframe*(10^(-rmc.PDSCH.Rho/20)),chSubframe,noiseest);

            % Transport block sizes
            outLen = ones(1,ncw)*trblksize(rmc.NSubframe+1);

            % Decode DownLink Shared Channel (DL-SCH)
            [decbits,blkcrc] = lteDLSCHDecode(rmc,rmc.PDSCH,outLen,...
                               rxEncodedBits);
                           
            % Skip empty transport blocks, this will be the case of
            % subframe 5, where no data is transmitted in RMC R.12
            if(outLen ~= 0) 
                trSizes(dataSubframeIndex, :) = outLen;
                crc(dataSubframeIndex, :) = blkcrc;
                dataSubframeIndex = dataSubframeIndex + 1;
            end
        end

        % Store resulting block CRC results and bit throughput vector
        % for analysis
        TotalBLKCRCcws(FrameNo,:) = crc(:);
        BitTput(FrameNo,:) = trSizes(:).*(1 - crc(:));
        rmc.NSubframe = 0;
    end
    
    % Record the block CRC and bit throughput for the total number of
    % frames simulated at a particular SNR
    TotalBLKCRC(resultIndex,:) = TotalBLKCRCcws(:); % CRC 
    BitThroughput(resultIndex,:) = BitTput(:); % successfully decoded bits
    resultIndex = resultIndex + 1;
    offsets = 0;
end

figure
plot(SNRIn,mean(BitThroughput,2),'-*')
title(['Throughput for ', num2str(NFrames) ' frame(s)'] )
xlabel('SNRdB'); ylabel('Throughput (kbps)');
grid on;
hold on;
plot(SNRIn,mean([trblksize(1:5) trblksize(7:end)])*0.7*...
    ones(1,numel(SNRIn)),'--rs');
legend('Simulation Result','70 Percent Throughput','Location',...
    'SouthEast')

figure
plot(SNRIn,100*(1-mean(TotalBLKCRC,2)),'-*')
title(['Throughput for ', num2str(NFrames) ' frame(s)'] )
xlabel('SNRdB'); ylabel('Throughput (%)');
grid on;
hold on;
plot(SNRIn,70*ones(1,numel(SNRIn)),'--rs');
legend('Simulation Result','70 Percent Throughput','Location',...
    'SouthEast')
ylim([0 100])

enb.PDSCH.TxScheme = 'CDD';         % Transmission scheme
enb.PDSCH.Rho = -6;                 % Downlink power allocation
enb.RC = 'R.14';                    % RMC number