Test = 1; % Test number (1 or 2)

if (Test==1)
    targetSNR = 2.6;
    SNRIn = [-1.5 1.0 3.5];
else
    targetSNR = -3.2;
    SNRIn = [-7.3 -4.8 2.3];
end

nFrames = 1;

% eNodeB configuration according to TS 36.101 Table 8.8.1.1-2
enb = struct;
enb.NDLRB = 50;
enb.DuplexMode = 'FDD';

% eNodeB configuration according to TS 36.101 Table 8.8.1.1-1
enb.CyclicPrefix = 'Normal';
enb.CellRefP = 2;
enb.CFI = 2 - (enb.NDLRB <= 10); % 2 symbols for PDCCH
enb.NCellID = 0;

% EPDCCH configuration according to TS 36.101 Table 8.8.1.1-1
epdcch = struct;
epdcch.EPDCCHTypeList = {'Distributed' 'Distributed'};
epdcch.EPDCCHPRBSetList = {[0; 7; 14; 21; 28; 35; 42; 49] [3; 17; 31; 45]};
epdcch.ControlChannelType = 'EPDCCH';
epdcch.DCIFormat = 'Format2A';

% EPDCCH configuration according to TS 36.101 Table 8.8.1.1-2
if (Test==1)
    epdcch.EPDCCHFormat = 1; % Aggregation level L = 4 ECCE
else
    epdcch.EPDCCHFormat = 3; % Aggregation level L = 16 ECCE
end

% Configure the EPDCCH scrambling identity and RNTI. These values are not
% mandated by the test configuration, but the parameters are required.
epdcch.EPDCCHNID = 0;
epdcch.RNTI = 1;

channel = struct;
channel.NRxAnts = 2;
channel.DelayProfile = 'EVA';
channel.MIMOCorrelation = 'Low';
if (Test==1)
    channel.DopplerFreq = 5;
else
    channel.DopplerFreq = 70;
end

perfectChannelEstimator = false;

if (~perfectChannelEstimator)
    cec.PilotAverage = 'UserDefined';   % Type of pilot symbol averaging
    cec.FreqWindow = 1;                 % Frequency window size in REs
    cec.TimeWindow = 2;                 % Time window size in REs
    cec.InterpType = 'linear';          % 2D interpolation type
    cec.InterpWindow = 'Centered';      % Interpolation window type
    cec.InterpWinSize = 1;              % Interpolation window size
    cec.Reference = 'EPDCCHDMRS';       % EPDCCH DM-RS used as reference
end

txSetIndex = 3-Test; % Test=1 -> 2nd set, Test=2 -> 1st set

% Configure EPDCCH set for transmission from the two EPDCCH sets according
% to 'txSetIndex'
epdcch.EPDCCHType = epdcch.EPDCCHTypeList{txSetIndex};
epdcch.EPDCCHPRBSet = epdcch.EPDCCHPRBSetList{txSetIndex};

% The 'mdsg' matrix will contain ones in the locations of missed downlink
% scheduling grants; the 'mcfi' matrix will contain ones in the locations
% of missed CFI decoding
mdsg = zeros(numel(SNRIn),10*nFrames);
mcfi = zeros(numel(SNRIn),10*nFrames);

epdcchports = 2;
portcounts = [enb.CellRefP epdcchports];
if (isfield(enb,'CSIRefP'))
    portcounts = [portcounts enb.CSIRefP];
end
ntxants = max(portcounts);

enb.NSubframe = 0;
[~,candidateInfo] = lteEPDCCHSpace(enb,epdcch);

fprintf('\n-- Parameter summary: --------------------------------------------\n');
fprintf('  Number of transmitter antennas: %d\n',ntxants);
fprintf('               Channel bandwidth: %gMHz\n',hNRBToBandwidth(enb.NDLRB));
fprintf('Number of OFDM symbols for PDCCH: %d symbols\n',enb.CFI + (enb.NDLRB <= 10));
fprintf('               Aggregation level: %d ECCE \n',candidateInfo.NECCEPerEPDCCH);
fprintf('                      DCI format: %s\n',epdcch.DCIFormat);
disp('------------------------------------------------------------------');

% For each SNR point:
for SNRIdx = 1:numel(SNRIn)

    % Set the random number generator to its default value
    rng('default');

    % Get current SNR value
    SNRdB = SNRIn(SNRIdx);

    fprintf('\nSimulating at %gdB SNR for a total of %d Frame(s)\n',SNRdB,nFrames);

    % Initialize timing offset vector
    offsets = 0;

    % For each subframe:
    for sf = 0:(nFrames*10)-1

        % Update frame and subframe numbers
        enb.NFrame = floor(sf/10);
        enb.NSubframe = mod(sf,10);

        % Configure EPDCCH candidate; the first candidate from the set of
        % search space candidates is used
        candidates = lteEPDCCHSpace(enb,epdcch);
        epdcch.EPDCCHECCE = candidates(1,:);

        % Create random DCI message bits from the configured DCI format
        dciInfo = lteDCIInfo(enb,epdcch);
        [dci,dciBits] = lteDCI(enb,epdcch,randi([0 1],dciInfo.(epdcch.DCIFormat),1));

        % Create EPDCCH indices; the 2nd output contains the EPDCCH bit
        % capacity EPDCCHG
        [epdcchIndices,epdcchInfo] = lteEPDCCHIndices(enb,epdcch);

        % Encode the DCI message bits to match the bit capacity EPDCCHG
        codedDciBits = lteDCIEncode(epdcch,dciBits,epdcchInfo.EPDCCHG);

        % Modulate the coded DCI message on the EPDCCH
        epdcchSymbols = lteEPDCCH(enb,epdcch,codedDciBits);

        % Create EPDCCH Demodulation Reference Signal (DM-RS)
        epdcchDmrsSymbols = lteEPDCCHDMRS(enb,epdcch);
        epdcchDmrsIndices = lteEPDCCHDMRSIndices(enb,epdcch);

        % Create empty subframe resource grid
        subframe = lteDLResourceGrid(enb,ntxants);

        % Concatenate EPDCCH symbols/indices with EPDCCH DM-RS
        % symbols/indices to facilitate beamforming
        allSymbols = [epdcchSymbols; epdcchDmrsSymbols];
        allIndices = [epdcchIndices; epdcchDmrsIndices];

        % Perform EPDCCH beamforming according to TS 36.101 Annex B.4.4 for
        % Distributed type and Annex B.4.5 for Localized type
        [K,L,P] = size(subframe);
        [resubs,~,portsubs] = ind2sub([K L 4],allIndices);
        rbsubs = floor((resubs-1)/12)+1;
        rbs = unique(rbsubs);
        ports = unique(portsubs.');
        if (strcmpi(epdcch.EPDCCHType,'Localized'))
            codebookIdx = randi([0 3],1);
            W = lteDLPrecode(1,ntxants,'SpatialMux',codebookIdx);
        end
        codebookIdxs = zeros(length(rbs),length(ports));
        for p=1:length(ports)
            thisport = (portsubs==ports(p));
            for r = 1:length(rbs)
                thisrb = (rbsubs==rbs(r));
                if (strcmpi(epdcch.EPDCCHType,'Distributed'))
                    unusedIdxs = setxor(0:3,codebookIdxs(r,1:p-1));
                    codebookIdx = unusedIdxs(randi(length(unusedIdxs),1));
                    W = lteDLPrecode(1,ntxants,'SpatialMux',codebookIdx);
                end
                codebookIdxs(r,p) = codebookIdx;
                bfSymbols = allSymbols(thisport & thisrb) * W;
                [~,bfIndices] = lteExtractResources(allIndices(thisport & thisrb),subframe);
                subframe(bfIndices) = subframe(bfIndices) + bfSymbols;
            end
        end

        % Transmit CFI on the PCFICH
        subframe(ltePCFICHIndices(enb)) = ltePCFICH(enb,lteCFI(enb));

        % Transmit Cell-Specific Reference Signals (CRS)
        subframe(lteCellRSIndices(enb)) = lteCellRS(enb);

        % Transmit PSS and SSS
        subframe(ltePSSIndices(enb)) = ltePSS(enb);
        subframe(lteSSSIndices(enb)) = lteSSS(enb);

        % Add OCNG
        subframe = addOCNG(enb,epdcch,subframe);

        % Perform OFDM modulation
        [txwaveform,info] = lteOFDMModulate(enb,subframe);

        % Set sampling rate of channel to that of the OFDM modulation
        channel.SamplingRate = info.SamplingRate;

        % Set channel time to current subframe
        channel.InitTime = sf/1000;

        % Reinitialize channel seed for each subframe to increase
        % variability
        channel.Seed = sf + 1;

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

        % Linear SNR
        SNR = 10^(SNRdB/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*ntxants*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
        if (enb.NSubframe==0)
            offset = lteDLFrameOffset(enb,rxwaveform);
            % Determine if frame synchronization was successful
            if (offset > 25)
                offset = offsets(end);
            else
                offsets = [offsets offset]; %#ok
            end
        end
        if (offset>0)
            rxwaveform = rxwaveform(1+offset:end,:);
        end

        % Perform OFDM demodulation
        rxsubframe = lteOFDMDemodulate(enb,rxwaveform);

        % Perform CRS-based (or perfect) channel estimation for PCFICH /
        % CFI reception
        if (perfectChannelEstimator)
            % Perform perfect channel estimation assuming the number of
            % transmit antennas is CellRefP (i.e. for CRS-based reception)
            hest = lteDLPerfectChannelEstimate(enb,channel,offset); %#ok<UNRCH>
            noiseGrid = lteOFDMDemodulate(enb,noise(1+offset:end ,:));
            nest = var(noiseGrid(:));
        else
            % Perform channel estimation using the CRS. The channel
            % estimator averaging parameters are adjusted to use a
            % rectangular averaging window suited to the CRS, rather than
            % the 1-by-2 (frequency-by-time) window configured in 'cec'
            % which is suited to the EPDCCH DM-RS symbols (which are
            % adjacent in time).
            cec_crs = cec;
            cec_crs.Reference = 'CellRS';
            cec_crs.FreqWindow = 31;
            cec_crs.TimeWindow = 15;
            [hest,nest] = lteDLChannelEstimate(enb,cec_crs,rxsubframe);
        end

        % Receive PCFICH / CFI
        pcfichIndices = ltePCFICHIndices(enb);
        [pcfichRx,pcfichHest] = lteExtractResources(pcfichIndices,rxsubframe,hest);
        rxCFI = lteCFIDecode(ltePCFICHDecode(enb,pcfichRx,pcfichHest,nest));

        % If the CFI was incorrectly decoded:
        if (rxCFI~=enb.CFI)

            % Record a misdetection of the CFI as a misdetection of the DCI
            % message (in 'mdsg'), in accordance with the test definition
            % in TS 36.101 Section 8.8. It is also recorded separately as a
            % misdetection of the CFI (in 'mcfi'). 'mcfi' can be used to
            % establish the probability of misdetection of the CFI on its
            % own. At low SNRs and for robust EPDCCH configurations (high
            % aggregation level), misdetection of the CFI may dominate the
            % overall misdetection of the DCI message.
            mdsg(SNRIdx,sf+1) = 1;
            mcfi(SNRIdx,sf+1) = 1;

        else % Proceed to EPDCCH / DCI decoding

            if (perfectChannelEstimator)
                if (ntxants~=enb.CellRefP) %#ok<UNRCH>
                    % Perform perfect channel estimation for 'ntxants'
                    % antennas (i.e. for EPDCCH reception) if this differs
                    % from CellRefP (used earlier for PCFICH/CFI reception)
                    hest = lteDLPerfectChannelEstimate(enb,channel,offset,ntxants);
                end
                % Apply EPDCCH beamforming to the perfect channel estimate
                hest = beamformPerfectChannelEstimate(hest,rbs,ports,codebookIdxs);
            end

            % Initialize array for received DCI message
            rxDCI = {};

            % Copy the structure 'epdcch' to 'epdcchRx' and remove fields
            % of 'epdcchRx' which are not required (they will be set during
            % decoding and search). This is not essential but provides
            % realism because |epdcchRx| then only contains system
            % information conveyed to the UE from the RRC
            epdcchRx = epdcch;
            epdcchRx = rmfield(epdcchRx,{'DCIFormat','EPDCCHFormat','EPDCCHECCE'});

            % For each EPDCCH set:
            for rxSetIndex = 1:numel(epdcchRx.EPDCCHTypeList)

                % Configure reception for current set
                epdcchRx.EPDCCHType = epdcchRx.EPDCCHTypeList{rxSetIndex};
                epdcchRx.EPDCCHPRBSet = epdcchRx.EPDCCHPRBSetList{rxSetIndex};

                if (~perfectChannelEstimator)
                    % Perform channel estimation using the EPDCCH DM-RS
                    [hest,nest] = lteDLChannelEstimate(enb,epdcchRx,cec,rxsubframe);
                end

                % Perform EPDCCH demodulation
                [bits,symbols] = lteEPDCCHDecode(enb,epdcchRx,rxsubframe,hest,nest);

                % Perform EPDCCH blind search and DCI message decoding
                rxDCI = [rxDCI lteEPDCCHSearch(enb,epdcchRx,bits)]; %#ok<AGROW>

            end

            % Record if the DCI message was missed; specifically if no DCI
            % messages were decoded or if none of the decoded messages
            % matched the transmitted message
            mdsg(SNRIdx,sf+1) = (isempty(rxDCI) || ~any(cellfun(@(x)isequal(x,dci),rxDCI)));

        end

    end

end

Pmdsg = 100*mean(mdsg,2);

plotResults(SNRIn,targetSNR,nFrames,Pmdsg);

function hest = beamformPerfectChannelEstimate(hestperfect,rbs,ports,codebookIdxs)

    [K,L,nrxants,ntxants] = size(hestperfect);
    hest = zeros([K L nrxants 4]);
    for r = 1:length(rbs)
        resubs = (rbs(r)-1)*12 + (1:12).';
        allportW = zeros(4,ntxants);
        activeportW = arrayfun(@(x)lteDLPrecode(1,ntxants,'SpatialMux',x),codebookIdxs(r,:),'UniformOutput',false);
        allportW(ports,:) = cat(1,activeportW{:});
        hest(resubs,:,:,:) = reshape(reshape(hestperfect(resubs,:,:,:),[12*L*nrxants ntxants])*allportW.',[12 L nrxants 4]);
    end

end

function subframe = addOCNG(enb,epdcch,subframe)

    % Add OP.7 OCNG (TS 36.101 Annex A.5.1.7). Consists of a single PDSCH
    % in all resource blocks not used by the EPDCCH.
    ocngPdsch.Modulation = 'QPSK';
    ocngPdsch.RNTI = 0;
    if (enb.CellRefP==1)
        ocngPdsch.TxScheme = 'Port0';
    else
        ocngPdsch.TxScheme = 'TxDiversity';
    end
    ocngPRBs = setdiff((0:enb.NDLRB-1).',epdcch.EPDCCHPRBSet);
    [ocngPdschInd,ocngPdschInfo] = ltePDSCHIndices(enb,ocngPdsch,ocngPRBs);
    ocngPdschData = randi([0 1],1,ocngPdschInfo.G);
    ocngPdschPower = -3.0; % rho + sigma = -3.0 + 0.0
    ocngPdschSym = ltePDSCH(enb,ocngPdsch,ocngPdschData)*(10^(ocngPdschPower/20));
    subframe(ocngPdschInd) = ocngPdschSym;

    % Add control region OCNG reference channel (TS 36.101 Annex A.5).
    % Includes PDCCH and PHICH locations but avoids PCFICH locations, to
    % allow the CFI to be transmitted
    pcfichRegInd = ltePCFICHIndices(enb,'reg');
    enbPdcch = enb;
    enbPdcch.PHICHDuration = 'Normal';
    enbPdcch.Ng = 'Sixth';
    exreg = pcfichRegInd(:,1);
    ocngPdcchInd = ltePDCCHIndices(enbPdcch,exreg);
    NRE = size(ocngPdcchInd,1);
    ocngPdcchData = randi([0 1],1,NRE*2);
    ocngPdcchPower = -3.0; % rho + sigma = -3.0 + 0.0
    ocngPdcchSym = ltePDCCH(enbPdcch,ocngPdcchData,NRE/4)*(10^(ocngPdcchPower/20));
    subframe(ocngPdcchInd) = ocngPdcchSym;

end

function plotResults(SNRIn,targetSNR,nFrames,Pmdsg)

    figure;
    plot(SNRIn,Pmdsg,'-*')
    if (nFrames==1)
        frameStr = '1 frame';
    else
        frameStr = [num2str(nFrames) ' frames'];
    end
    title(sprintf('Probability of missed downlink scheduling grant\nfor %s',frameStr));
    xlabel('SNR (dB)'); ylabel('Pm-dsg (%)');
    grid on;
    hold on;
    plot(targetSNR,1,'--r^','MarkerFaceColor','r')
    plot(SNRIn, ones(1,numel(SNRIn)),'--r');
    legend('Simulation Result',sprintf('Target: 1%% Pm-dsg @ %0.1fdB SNR',targetSNR),'Location','NorthEast');
    xlim([min(SNRIn)-0.5 max(SNRIn)+0.5]);
    ylim([fix(min(Pmdsg)-2) fix(max(Pmdsg)+2)]);

end