DL_model.m

%% uplink model for NB-IoT NTN

function [] = DL_model(SAT_EIRP, G_T, SAT_altitude)
 
    satelliteParamsSource = "Custom";
    %% Satellite parameters
    % When you set satelliteParamsSource to Custom, provide these fields.
    satellite = struct;
    satellite.EIRPDensity = SAT_EIRP;   % dBW/MHz
    satellite.RxGByT = G_T;                 % dB/K
    satellite.Altitude = SAT_altitude;      % 60000m = 600km
    
    %% Set User Equipment Payload Characteristics
    % By default, this example uses the UE power class PC3 with a noise figure of 7 dB
    ue = struct;
    ue.TxPower = 23;               % dBm
    ue.TxGain = 0;                 % dBi
    ue.TxCableLoss = 0;            % dB
    ue.RxNoiseFigure = 7;          % dB
    ue.RxGain = 0;                 % dBi
    ue.RxAntennaTemperature = 290; % K
    ue.RxAmbientTemperature = 290; % K
    
    % Set Link Characteristics
    % Set the link characteristics by defining the link direction, elevation angle, bandwidth, frequency, and atmospheric losses.
    % The transmission bandwidths for NB-IoT are:
    % Downlink — 180 kHz
    % Uplink — 3.75 kHz, 15 kHz, 45 kHz, 90 kHz, or 180 kHz
    
    % Enable or disable the ITU-R P.618 propagation losses
    
    useP618PropagationLosses = false;           % 0 (false), 1 (true)
    
    %% Set the link characteristics
    link = struct;
    link.Direction = "downlink";    % "uplink", "downlink"
    link.ElevationAngle = [5 15 30 60 90];         % degrees (Vector of elevation angles)
    link.Frequency = 2e9;                     % Hz
    link.Bandwidth = 180e3;                   % Hz
    link.ShadowMargin = 3;                    % dB
    link.AdditionalLosses = 0;                % dB
    link.PolarizationLoss = 0;                % dB
    
    if useP618PropagationLosses == 0
        % Set these fields when you set useP618PropagationLosses to false.
        link.ScintillationLosses = 2.2;           % dB
        link.AtmosphericLosses = 0.07;             % dB
    else
        % Set these fields when you set useP618PropagationLosses to true. In
        % this case, the example uses digital maps to calculate the
        % attenuation.
        link.P618Configuration = p618Config;
        link.P618Configuration.Latitude = 51.5;                   % degrees
        link.P618Configuration.Longitude = -0.14;                 % degrees
        link.P618Configuration.GasAnnualExceedance = 1;
        link.P618Configuration.CloudAnnualExceedance = 1;
        link.P618Configuration.ScintillationAnnualExceedance = 1;
        link.P618Configuration.TotalAnnualExceedance = 1;
        link.P618Configuration.PolarizationTiltAngle = 0;         % degrees
        link.P618Configuration.AntennaDiameter = 1;               % meters
        link.P618Configuration.AntennaEfficiency = 0.5;
    end
    
    %% Set NB-IoT Data Channel Parameters
    % To calculate the reference carrier-to-noise ratio (CNR) for an NB-IoT data channel transmission
    % set these NB-IoT parameters
    % - Modulation scheme
    % - Transport block size (Nbits)
    % - Number of symbols used for transmission (Nsymbols)
    % - Number of repetitions (NRep)
    % - Number of subframes in downlink (NSF)
    % - Number of resource units in uplink (NRU)
    
    modulation = "QPSK"; % Modulation scheme
    nBits = 120;                       % TBS, Number of useful bits (transport block size)
    nSymbols = 160;                    % Number of symbols used for transmission
    nRep = 2;                          % Number of repetitions
    nSF = 8;                           % Number of subframes (used in downlink)
    nRU = 1;                           % Number of resource units (used in uplink)
    
    %% Set these NB-IoT waveform parameters.
    % - Number of used data subcarriers (NDSC)
    % - Number of fast Fourier transform (FFT) bins (NFFT)
    % - Data symbol duration (Td)
    % - Cyclic prefix duration (TCP)
    % - Oversampling factor (OSR)
    
    nDSC = 72;                         % Number of data subcarriers (6 resource blocks)
    nFFT = 128;                        % FFT length
    tsamp = 1/1.92e6;                  % Sample time (in s)
    td = nFFT*tsamp;                   % Data symbol duration (in s)
    nCP = 9;                           % Number of cyclic prefix samples
    tCP = nCP*tsamp;                   % Cyclic prefix duration (in s)
    osr = 1;                           % Oversampling factor
    
    % In the presence of additive white Gaussian noise without any channel coding
    % the modulation schemes return these results:
    % - BPSK achieves a bit error rate of 1e-6 at an Eb/No of 10.5dB
    % - QPSK achieves a bit error rate of 1e-6 at an Eb/No of 10.5dB
    % - 16-QAM achieves a bit error rate of 1e-6 at an Eb/No of 14.4dB
    
    %% Reference Eb/No in dB
    if modulation == "16-QAM"
        % Modulation scheme is 16-QAM
        ebnoRef = 14.4;
    else 
        % Modulation scheme is either BPSK or QPSK
        ebnoRef = 10.5;
    end
    
    %% Compute Reference CNR
    % To compute the reference CNR in dB, use this equation:
    % (C/N)ref = (Eb/No)ref + 10*log(m*Reff) + 10*log(Ndsc/Nfft) + 10*log(Td/(Td+Tcp)) - 10*log(OSR)
    
    % (Eb/No)ref is the value in dB for the desired bit error rate
    %　m is the modulation order
    % Reff is the effective code rate
    
    %% The effective code rate for NB-IoT downlink:
    % Reff = (nBits + nCRC)/(nSF*nSymbols*m*nRep)
    % - Nbits is the number of useful bits to transmit (transport block size)
    % - NCRC is the number of bits used for check code (24 for NB-IoT)
    % - NSF is the number of allocated subframes
    % - Nsymbols is the number of symbols
    % - Nrep is the number of repetitions
    
    %% The effective code rate for NB-IoT uplink:
    % Reff = (nBits + nCRC)/(nRU*nSymbols*m*nRep)
    % NRU is the number of allocated resource units
    % For 1-tone transmission, Nsymbols is 96
    % for 3-tone transmission, Nsymbols is 144
    
    % Modulation order
    m = 1;
    if modulation == "QPSK"
        m = 2;
    elseif modulation == "16-QAM"
        m = 4;
    end
    
    % Number of CRC bits
    nCRC = 24;
    
    % Calculate the effective code rate and the reference CNR
    if link.Direction == "downlink"
        Reff = (nBits + nCRC)/(nSF*nSymbols*m*nRep);
    else % "uplink"
        Reff = (nBits + nCRC)/(nRU*nSymbols*m*nRep);
    end
    cnrRef = ebnoRef + 10*log10(m*Reff) + 10*log10(nDSC/nFFT) + ...
            10*log10(td/(td + tCP)) - 10*log10(osr);
    
    disp("Reference CNR: " + cnrRef + " dB")
    
    %% Compute Link Budget
    % EIRP is the effective isotropic radiated power in dBW.
    % G/T is the antenna-gain-to-noise temperature in dB/K.
    % k is the Boltzmann constant with the value of -228.6 dBW/K/Hz.
    % PL_FS is the free space path loss (FSPL) in dB.
    % PL_A is the atmospheric path loss due to gases and shadowing margin in dB.
    % PL_SM is the shadowing margin in dB.
    % PL_SL is the scintillation loss in dB.
    % PL_AD is the additional loss in dB.
    % B is the channel bandwidth in dBHz.
    
    % To calculate antenna-gain-to-noise-temperature
    % GR is the receive antenna gain in dBi.
    % Nf is the noise figure in dB.
    % T0 is the ambient temperature in degrees Kelvin.
    % Ta is the antenna temperature in degrees Kelvin.
    
    % The receive antenna gain depends on the type of antenna used, and accounts for polarization loss.
    % EIRP = PT - LC + GT
    % PT is the transmit antenna power in dBW.
    % LC is the cable loss in dB.
    % GT is the transmit antenna gain in dBi.
    
    % Get the satellite parameters based on satelliteParamsSource
    if satelliteParamsSource ~= "Custom"
        satellite = getSatelliteParams(satelliteParamsSource);
    end
    
    % Calculate EIRP when the satellite is a transmitter (dB)
    % To get a value in decibels from density: 
    % Value in dB = Value in dB/MHz + 10*log10[BW in MHz]
    satellite.EIRP = satellite.EIRPDensity + 10*log10(link.Bandwidth/1e6);
    
    % Calculate EIRP when the UE is a transmitter (dB)
    % To get a value in decibels from dBm: dB = dBm - 30
    ue.EIRP = (ue.TxPower-30) + ue.TxGain - ue.TxCableLoss;
    
    % Calculate the gain-to-noise temperature or the figure of merit for the UE
    % as a receiver.
    ue.RxGByT = ue.RxGain - ue.RxNoiseFigure ...
        - 10*log10(ue.RxAmbientTemperature + ...
        (ue.RxAntennaTemperature-ue.RxAmbientTemperature)*10^(-0.1*ue.RxNoiseFigure));
    
    % Set the transmitter and the receiver based on the link direction
    if link.Direction == "uplink"
        tx = ue;
        rx = satellite;
    else
        tx = satellite;
        rx = ue;
    end
    
    % Range of elevation angles
    elevAngles = link.ElevationAngle(:);
    numElevAngles = numel(elevAngles);
    
    % Calculate the distance from the satellite to the ground station for all
    % elevation angles
    re = physconst("earthradius");
    c = physconst("lightspeed");
    h = satellite.Altitude;
    d = -re*sind(elevAngles) + sqrt((re*sind(elevAngles)).^2 + h*h + 2*re*h);
    
    % Get the total atmospheric losses for each elevation angle
    totalAtmosphericLoss = zeros(numElevAngles,1);
    if useP618PropagationLosses == 1
        % Download and extract the digital maps, if not available on the path
        maps = exist("maps.mat","file");
        p836 = exist("p836.mat","file");
        p837 = exist("p837.mat","file");
        p840 = exist("p840.mat","file");
        matFiles = [maps p836 p837 p840];
        if ~all(matFiles)
            if ~exist("ITURDigitalMaps.tar.gz","file")
                url = "https://www.mathworks.com/supportfiles/spc/P618/ITURDigitalMaps.tar.gz";
                websave("ITURDigitalMaps.tar.gz",url);
                untar("ITURDigitalMaps.tar.gz")
            else
                untar("ITURDigitalMaps.tar.gz")
            end
        end
        link.P618Configuration.Frequency = link.Frequency;
        elevAnglesToConsider = elevAngles;
        if any(elevAngles < 5)
            warning("The prediction method for scintillation losses is valid for elevation " + ... 
                "angle greater than 5 degree. For elevation angle less than 5 degree, the " + ... 
                "nearest valid value of 5 degree will be used in the computation.")
            elevAnglesToConsider(elevAngles < 5) = 5;
        end
        for index = 1:numElevAngles
            link.P618Configuration.ElevationAngle = elevAnglesToConsider(index);
            pl = p618PropagationLosses(link.P618Configuration);
            totalAtmosphericLoss(index) = pl.At;
        end
    else
        totalAtmosphericLoss(:) = link.AtmosphericLosses + link.ScintillationLosses;
    end
    
    % Define the configuration parameters
    config = satelliteCNRConfig;
    config.TransmitterPower = tx.EIRP;
    config.TransmitterAntennaGain = 0;
    config.Frequency = link.Frequency/1e9;          % GHz
    config.GainToNoiseTemperatureRatio = rx.RxGByT;
    config.Bandwidth = link.Bandwidth/1e6;          % MHz
    
    % Compute the CNR and the free space path loss for each elevation angle
    cnr = zeros(numElevAngles,1);
    pathLoss = cnr;
    for index = 1:numElevAngles
        config.Distance = d(index)/1e3;                                         % km
        config.MiscellaneousLoss = totalAtmosphericLoss(index) + ...
             link.PolarizationLoss + link.ShadowMargin + link.AdditionalLosses;
        [cnr(index),cnrInfo] = satelliteCNR(config);
        pathLoss(index) = cnrInfo.FSPL;
    end
    
    disp("Downlink ==========")
    % Report the results in a table
    table(elevAngles,cnr,pathLoss,VariableNames=["Elevation Angle (degrees)", "CNR (dB)", "FSPL (dB)"])
    
    % Plot CNR as a function of the elevation angle
    plot(elevAngles,cnr,"*")
    title("CNR As a Function of Elevation Angle")
    xlabel("Elevation Angle (degrees)")
    ylabel("CNR (dB)")
    ylim([min(cnr)-0.2 max(cnr)+0.2])
    xlim([min(elevAngles)-1 max(elevAngles)+1])
    grid on
    
    %% Compute Link Margin and NB-IoT Repetitions
    % The link margin (in dB) is the difference between the calculated CNR and the reference CNR. 
    % The link is closed when the link margin is zero or positive.
    % To calculate the link margin, use this equation: 
    % LinkMargin = C/N-(C/N)ref
    
    % Compute link margin
    linkMargin = cnr - cnrRef;
    
    %% Minimum number of additional repetitions required for link closure
    minRepetitions = 10.^(-linkMargin./10);
    idx = linkMargin >= 0;
    additionalRepetitions = minRepetitions;
    additionalRepetitions(idx) = 0;
    % When the link margin is negative, improve the CNR by adding repetitions.
    % Calculate the required number of additional repetitions.
    additionalRepetitions(~idx) = ceil(nRep*(additionalRepetitions(~idx)-1));
    
    % Report the results in a table with this format.
    % Elevation Angle | Link Margin | Min. Additional Repetitions Required
    table(elevAngles,linkMargin,additionalRepetitions, ...
        VariableNames=["Elevation Angle (degrees)", ...
        "Link Margin (dB)", "NRep_Add"])
    %% Local Function
    % getSatelliteParams — Gets the satellite parameters defined for the selected satelliteParamsSource
    function satellite = getSatelliteParams(satelliteParamsSource)
    
        if satelliteParamsSource == "Set 1 GEO"
            % Table 6.2-4, TR 36.763
            satellite = struct;
            satellite.EIRPDensity = 59;    % dBW/MHz
            satellite.RxGByT = 19;         % dB/K
            satellite.Altitude = 35786e3;  % m
        elseif satelliteParamsSource == "Set 1 LEO-1200"
            % Table 6.2-4, TR 36.763
            satellite = struct;
            satellite.EIRPDensity = 40;    % dBW/MHz
            satellite.RxGByT = 1.1;        % dB/K
            satellite.Altitude = 1200e3;   % m
        elseif satelliteParamsSource == "Set 1 LEO-600"
            % Table 6.2-4, TR 36.763
            satellite = struct;
            satellite.EIRPDensity = 34;    % dBW/MHz
            satellite.RxGByT = 1.1;        % dB/K
            satellite.Altitude = 600e3;    % m
        elseif satelliteParamsSource == "Set 2 GEO"
            % Table 6.2-5, TR 36.763
            satellite = struct;
            satellite.EIRPDensity = 53.5;  % dBW/MHz
            satellite.RxGByT = 14;         % dB/K
            satellite.Altitude = 35786e3;  % m
        elseif satelliteParamsSource == "Set 2 LEO-1200"
            % Table 6.2-5, TR 36.763
            satellite = struct;
            satellite.EIRPDensity = 34;    % dBW/MHz
            satellite.RxGByT = -4.9;       % dB/K
            satellite.Altitude = 1200e3;   % m
        elseif satelliteParamsSource == "Set 2 LEO-600"
            % Table 6.2-5, TR 36.763
            satellite = struct;
            satellite.EIRPDensity = 28;    % dBW/MHz
            satellite.RxGByT = -4.9;       % dB/K
            satellite.Altitude = 600e3;    % m
        elseif satelliteParamsSource == "Set 3 GEO"
            % Table 6.2-6, TR 36.763
            satellite = struct;
            satellite.EIRPDensity = 59.8;  % dBW/MHz
            satellite.RxGByT = 16.7;       % dB/K
            satellite.Altitude = 35786e3;  % m
        elseif satelliteParamsSource == "Set 3 LEO-1200"
            % Table 6.2-6, TR 36.763
            satellite = struct;
            satellite.EIRPDensity = 33.7;  % dBW/MHz
            satellite.RxGByT = -12.8;      % dB/K
            satellite.Altitude = 1200e3;   % m
        elseif satelliteParamsSource == "Set 3 LEO-600"
            % Table 6.2-6, TR 36.763
            satellite = struct;
            satellite.EIRPDensity = 28.3;  % dBW/MHz
            satellite.RxGByT = -12.8;      % dB/K
            satellite.Altitude = 600e3;    % m
        elseif satelliteParamsSource == "Set 4 LEO-600"
            % Table 6.2-7, TR 36.763
            satellite = struct;
            satellite.EIRPDensity = 21.45; % dBW/MHz
            satellite.RxGByT = -18.6;      % dB/K
            satellite.Altitude = 600e3;    % m
        else % "Set 5 MEO-10000"
            % Table 6.2-8, TR 36.763
            satellite = struct;
            satellite.EIRPDensity = 45.4;  % dBW/MHz
            satellite.RxGByT = 3.8;        % dB/K
            satellite.Altitude = 10000e3;  % m
        end
    end 
end