NSubframes = 10;                            % Number of subframes
SNRIn = [-10.2 -8.2 -6.2 -4.2 -2.2];        % SNR range in dB
NTxAnts = 1;                                % Number of transmit antennas

ue = struct;                                % UE config structure
ue.NULRB = 6;                               % 6 resource blocks (1.4 MHz)
ue.CyclicPrefixUL = 'Extended';             % Extended cyclic prefix
ue.Hopping = 'Off';                         % No frequency hopping
ue.NCellID =  9;
ue.Shortened = 0;                           % No SRS transmission
ue.NTxAnts = NTxAnts;

% Hybrid ARQ indicator bit set to one. Only one bit is required for PUCCH
% 1a
ACK = 1;
pucch = struct;                         % PUCCH config structure
% Set the size of resources allocated to PUCCH format 2. This affects the
% location of PUCCH 1 transmission
pucch.ResourceSize = 0;
pucch.DeltaShift = 1;                       % Delta shift PUCCH parameter
% Number of cyclic shifts used for PUCCH format 1 in resource blocks with a
% mixture of formats 1 and 2. This is the N1cs parameter
pucch.CyclicShifts = 0;
% Vector of PUCCH resource indices, one per transmission antenna. This is
% the n2pucch parameter
pucch.ResourceIdx = 0:ue.NTxAnts-1;

channel = struct;                       % Channel config structure
channel.NRxAnts = 2;                    % Number of receive antennas
channel.DelayProfile = 'ETU';           % Channel delay profile
channel.DopplerFreq = 70.0;             % Doppler frequency in Hz
channel.MIMOCorrelation = 'Low';        % Low MIMO correlation
channel.NTerms = 16;                    % Oscillators used in fading model
channel.ModelType = 'GMEDS';            % Rayleigh fading model type
channel.Seed = 13;                      % Channel seed
channel.InitPhase = 'Random';           % Random initial phases
channel.NormalizePathGains = 'On';      % Normalize delay profile power
channel.NormalizeTxAnts = 'On';         % Normalize for transmit antennas

% SC-FDMA modulation information: required to get the sampling rate
scfdmaInfo = lteSCFDMAInfo(ue);
channel.SamplingRate = scfdmaInfo.SamplingRate;   % Channel sampling rate

cec = struct;                     % Channel estimation config structure
cec.PilotAverage = 'UserDefined'; % Type of pilot averaging
cec.FreqWindow = 12;              % Frequency averaging window in REs (special mode)
cec.TimeWindow = 1;               % Time averaging window in REs (Special mode)
cec.InterpType = 'cubic';         % Cubic interpolation

% Preallocate memory for probability of detection vector
PMISS = zeros(size(SNRIn));

for nSNR = 1:length(SNRIn)

    % Missed or incorrect ACK detection counter
    missCount = 0;

    % Noise configuration
    SNR = 10^(SNRIn(nSNR)/20);              % Convert dB to linear
    % The noise added before SC-FDMA demodulation will be amplified by the
    % IFFT. The amplification is the square root of the size of the IFFT.
    % To achieve the desired SNR after demodulation the noise power is
    % normalized by this value. In addition, because real and imaginary
    % parts of the noise are created separately before being combined into
    % complex additive white Gaussian noise, the noise amplitude must be
    % scaled by 1/sqrt(2*ue.NTxAnts) so the generated noise power is 1
    N = 1/(SNR*sqrt(double(scfdmaInfo.Nfft)))/sqrt(2.0*ue.NTxAnts);
    % Set the type of random number generator and its seed to the default
    % value
    rng('default')

    % Loop for subframes
    offsetused = 0;
    for nsf = 1:NSubframes

        % Create resource grid
        ue.NSubframe = mod(nsf-1,10);
        reGrid = lteULResourceGrid(ue);

        % Create PUCCH 1 and its DRS
        pucch1Sym = ltePUCCH1(ue,pucch,ACK);
        pucch1DRSSym = ltePUCCH1DRS(ue,pucch);

        % Generate indices for PUCCH 1 and its DRS
        pucch1Indices = ltePUCCH1Indices(ue,pucch);
        pucch1DRSIndices = ltePUCCH1DRSIndices(ue,pucch);

        % Map PUCCH 1 and PUCCH 1 DRS to the resource grid
        reGrid(pucch1Indices) = pucch1Sym;
        reGrid(pucch1DRSIndices) = pucch1DRSSym;

        % SC-FDMA modulation
        txwave = lteSCFDMAModulate(ue,reGrid);

        % Channel state information: set the init time to the correct value
        % to guarantee continuity of the fading waveform
        channel.InitTime = (nsf-1)/1000;

        % Channel modeling
        % The additional 25 samples added to the end of the waveform are to
        % cover the range of delays expected from the channel modeling (a
        % combination of implementation delay and channel delay spread)
        rxwave = lteFadingChannel(channel,[txwave; zeros(25,ue.NTxAnts)]);

        % Add noise at receiver
        noise = N * complex(randn(size(rxwave)),randn(size(rxwave)));
        rxwave = rxwave + noise;

        % Receiver

        % Synchronization
        % An offset within the range of delays expected from the channel
        % modeling (a combination of implementation delay and channel
        % delay spread) indicates success
        offset = lteULFrameOffsetPUCCH1(ue,pucch,rxwave);
        if (offset < 25)
            offsetused = offset;
        end

        % SC-FDMA demodulation
        rxgrid = lteSCFDMADemodulate(ue,rxwave(1+offsetused:end,:));

        % Channel estimation
        [H,n0] = lteULChannelEstimatePUCCH1(ue,pucch,cec,rxgrid);

        % Extract REs corresponding to the PUCCH 1 from the given subframe
        % across all receive antennas and channel estimates
        [pucch1Rx,pucch1H] = lteExtractResources(pucch1Indices,rxgrid,H);

        % MMSE equalization
        eqgrid = lteULResourceGrid(ue);
        eqgrid(pucch1Indices) = lteEqualizeMMSE(pucch1Rx,pucch1H,n0);

        % PUCCH 1 demodulation/decoding
        rxACK = ltePUCCH1Decode(ue,pucch,length(ACK),eqgrid(pucch1Indices));

        % Detect missed (empty rxACK) or incorrect HARQ-ACK (compare
        % against transmitted ACK
        if (isempty(rxACK) || any(rxACK ~= ACK))
            missCount = missCount + 1;
        end

    end

    PMISS(nSNR) = (missCount/NSubframes);

end

plot(SNRIn,PMISS,'b-o','LineWidth',2,'MarkerSize',7);
hold on;
plot(-4.2,0.01,'rx','LineWidth',2,'MarkerSize',7);
xlabel('SNR (dB)');
ylabel('Probability of missed ACK detection');
title('ACK missed detection test (TS36.104 Section 8.3.2.1)');
axis([SNRIn(1)-0.1 SNRIn(end)+0.1 -0.05 .35]);
legend('simulated performance','target');