%% PUSCH HARQ-ACK Detection Conformance Test
% This example shows how to measure the false detection probability and
% missed detection probability of Hybrid Automatic Repeat Request ACK
% (HARQ-ACK) multiplexed on the Physical Uplink Shared Channel (PUSCH)
% using the LTE Toolbox(TM) under conformance test conditions as defined in
% TS36.104, Section 8.2.4.1 [ <#9 1> ].

% Copyright 2009-2018 The MathWorks, Inc.

%% Introduction
% Uplink Control Information (UCI) can be carried on the Physical Uplink
% Shared Channel (PUSCH) and may contain Hybrid Automatic Repeat Request
% ACK (HARQ-ACK) information. This example uses the LTE Toolbox to perform
% the "HARQ-ACK Multiplexed on PUSCH" conformance test specified in
% TS36.104, Section 8.2.4.1 [ <#9 1> ].
%
% Two performance requirements are defined for HARQ-ACK multiplexed on the
% PUSCH:
% 
% * _ACK false detection probability_ is the probability that ACK is
% detected when data is only sent on symbols where HARQ-ACK information can
% be allocated.
% * _The ACK missed detection probability_ is the conditional probability
% of not detecting an ACK when it was sent on PUSCH resources.
%
% The conditions specified for this test are two receive antennas, normal
% cyclic prefix, ETU70 channel, using Fixed Reference Channel (FRC) A4-3.
% The target defined in TS 36.104, Section 8.2.4.1 [ <#9 1> ] for 1.4 MHz
% bandwidth (6RB) is a false detection probability and missed detection
% probability of 1% at 13.8 dB.
%
% The test is run on a subframe by subframe basis at each SNR test point.
% For each subframe an uplink waveform is generated with and without ACK
% information, passed through a fading channel and the HARQ-ACK decoded.
% False detections and missed detections are recorded and the probability
% of error displayed for the range of SNRs tested.

%% Simulation Configuration
% The example is executed for a simulation length of 1 frame at a number of
% SNRs including the required 13.8 dB at which the test requirements for
% false and missed detection rates (1% in each case) are to be achieved as
% per TS36.104, Table 8.2.4.1-1 [ <#9 1> ]. A large number of
% |numSubframes| should be used to produce meaningful throughput results.
% |SNRdB| can be an array of values or a scalar.

numSubframes = 10;  % Number of frames to simulate at each SNR
SNRdB = [4.8 7.8 10.8 13.8 16.8]; % SNR points to simulate

%% UE Configuration
% To configure the transmitter, a few desired parameter fields are set in
% the structure |frc|, which is then passed to <matlab:doc('lteRMCUL')
% lteRMCUL> which will set all the other required parameter fields. At this
% point, the structure |frc| represents the complete configuration of the
% required transmission.

frc.TotSubframes = 1;   % Total number of subframes to generate
frc.NCellID = 10;       % Cell identity
frc.RC = 'A4-3';        % FRC number

% Populate FRC configuration structure with default values for A4-3
frc = lteRMCUL(frc);

%% Propagation Channel Model Configuration
% The propagation channel is configured via the structure |chcfg|, with
% settings per the test requirements. Note that the sampling rate of the
% propagation is determined from the sampling rate of the transmitted
% waveform, which can be established using <matlab:doc('lteSCFDMAInfo')
% lteSCFDMAInfo>.

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

% Set channel model sampling rate
info = lteSCFDMAInfo(frc);
chcfg.SamplingRate = info.SamplingRate;

%% Channel Estimator Configuration
% The channel estimator is configured using the structure |cec|. An ETU
% delay profile causes the channel to change quickly over time. Therefore a
% small pilot averaging frequency window of 9 Resource Elements (REs) is
% used. The Demodulation Reference Signal (DRS) is contained in only one
% symbol per slot, therefore a time averaging window of 1 RE is used. This
% will not include any pilots from an adjacent slot when averaging.

cec.PilotAverage = 'UserDefined'; % Type of pilot averaging
cec.FreqWindow = 9;               % Frequency averaging windows in REs
cec.TimeWindow = 1;               % Time averaging windows in REs
cec.InterpType = 'cubic';         % Interpolation type
cec.Reference = 'Antennas';       % Reference for channel estimation 

%% Loop for SNR Values
% A loop is used to run the simulation for a set of SNR points, given by
% the vector |SNRdB|.
%
% 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
% ($N_{\mathrm{FFT}}$). To ensure that the power of the additive noise is
% normalized after demodulation to achieve the desired SNR, the desired
% noise power is divided by $N_{\mathrm{FFT}}$. In addition, as real and
% imaginary parts of the noise are created separately before being combined
% into complex Additive White Gaussian Noise (AWGN), the noise amplitude
% must be scaled by $1/\sqrt{2}$ so the generated noise power is 1.
%
% At each SNR test point the probability of false and missed detection is
% calculated. These probabilities are calculated on a subframe by subframe
% basis using the following steps:
%
% * _Create Transmit Waveform:_ The Uplink Reference Measurement Channel
% (RMC) Tool <matlab:doc('lteRMCULTool') lteRMCULTool> is used to generate
% an uplink waveform containing random transport data and a HARQ-ACK bit
% in every odd subframe.
%
% * _Noisy Channel Modeling:_ The waveform is passed through a fading
% channel and AWGN added.
%
% * _Perform Synchronization and SC-FDMA Demodulation:_ The received
% symbols are synchronized to account for a combination of implementation
% delay and channel delay spread. The symbols are then SC-FDMA demodulated.
%
% * _Perform Channel and Noise Power Spectral Density Estimation and
% Equalization:_ The channel and noise power spectral density are estimated
% and the received PUSCH symbols equalized.
%
% * _Decode the PUSCH:_ The UpLink Shared Channel (UL-SCH) coding is
% determined and used to decode the PUSCH to recover the interleaved
% UL-SCH.
%
% * _Recover UCI ACK Bit:_ Deinterleave the PUSCH and decode the HARQ-ACK
% information to recover the UCI HARQ-ACK bit and detect HARQ-ACK false or
% missed detection based on whether an HARQ-ACK bit was encoded in the
% transmitted subframe or not.

% Initialize variables used in the simulation and analysis
pFalse = zeros(size(SNRdB));  % Probability of false detection at each SNR
pMissed = zeros(size(SNRdB)); % Probability of missed detection at each SNR

for nSNR = 1:length(SNRdB)

    % Initialize the random number generator stream
    rng('default'); 
       
    % Extract SNR to test
    SNR = 10^(SNRdB(nSNR)/20);
    
    % Scale noise to ensure the desired SNR after SC-FDMA demodulation
    N = 1/(SNR*sqrt(double(info.Nfft)))/sqrt(2.0);    

    offsetUsed = 0;
    falseCount = 0;
    falseTests = 0;
    missCount = 0;
    missTests = 0;

    for subframeNo = 0:(numSubframes-1)

        % Updating subframe number
        frc.NSubframe = mod(subframeNo, 10);

        % Transmit ACK on every odd subframe
        if (mod(subframeNo, 2)==0)
            ACK = [];
            falseTests = falseTests + 1;
        else
            ACK = 1;
            missTests = missTests + 1;
        end            

        % Create random data to transmit
        trblklen = frc.PUSCH.TrBlkSizes(frc.NSubframe+1);
        trblk = randi([0 1], trblklen, 1);                                               

        % Transmit a waveform with an additional 25 samples to cover the
        % range of delays expected from the channel modeling (a
        % combination of implementation delay and channel delay spread)
        txWaveform = [lteRMCULTool(frc, trblk, [], [], ACK); zeros(25, 1)];        

        % Pass waveform through fading channel model
        chcfg.InitTime = subframeNo/1000;
        rxWaveform = lteFadingChannel(chcfg, txWaveform);

        % Add noise
        noise = N*complex(randn(size(rxWaveform)), ...
            randn(size(rxWaveform)));
        rxWaveform = rxWaveform + noise;

        % Synchronization
        offset = lteULFrameOffset(frc, frc.PUSCH,rxWaveform);
        if (offset < 25)
            offsetUsed = offset;
        end

        % SC-FDMA demodulation
        rxSubframe = lteSCFDMADemodulate(frc, ...
            rxWaveform(1+offsetUsed:end, :));

        % Channel Estimation
        [estChannelGrid, noiseEst] = lteULChannelEstimate( ...
            frc, frc.PUSCH, cec, rxSubframe);

        % PUSCH indices for given subframe
        puschIndices = ltePUSCHIndices(frc,frc.PUSCH);

        % Minimum Mean Squared Error (MMSE) equalization
        rxSymbols = lteEqualizeMMSE(rxSubframe(puschIndices), ...
                        estChannelGrid(puschIndices), noiseEst);

        % Obtain UL-SCH coding information for current transport block and
        % ACK, and concatenate this information with the PUSCH / UL-SCH
        % configuration
        frc.PUSCH = lteULSCHInfo(frc, frc.PUSCH, trblklen, ...
            0, 0, 1, 'chsconcat');

        % Perform deprecoding, demodulation and descrambling on the
        % received data
        rxEncodedBits = ltePUSCHDecode(frc, frc.PUSCH, rxSymbols); 

        % UL-SCH channel deinterleaving
        [deinterleavedData, ccqi, cri, cack] = ...
            lteULSCHDeinterleave(frc, frc.PUSCH, rxEncodedBits);

        % HARQ-ACK decoding
        rxACK = lteACKDecode(frc.PUSCH, cack);
        
        % Detect false or missed HARQ-ACK
        if (isempty(ACK) && ~isempty(rxACK))
            falseCount = falseCount + 1;
        end            
        if (~isempty(ACK) && ~isequal(ACK,rxACK)) 
            missCount = missCount + 1;
        end       
          
    end     

    % Calculate false or missed HARQ-ACK probability
    pFalse(nSNR) = falseCount/falseTests;
    pMissed(nSNR) = missCount/missTests;

end

%% Analysis
% Plot the simulated results against the target performance as stipulated
% in the standard using <matlab:edit('hHARQACKResults.m')
% hHARQACKResults.m>.

hHARQACKResults(SNRdB, pFalse, pMissed);

%% Appendix
% This example uses the helper functions:
%
% * <matlab:edit('hHARQACKResults.m') hHARQACKResults.m>

%% Selected Bibliography
% # 3GPP TS 36.104 "Base Station (BS) radio transmission and reception"
