Main Content

NR NTN PDSCH Throughput

This example shows how to measure the physical downlink shared channel (PDSCH) throughput of a 5G New Radio (NR) link in a non-terrestrial network (NTN) channel, as defined by the 3GPP NR standard. The example implements the PDSCH and downlink shared channel (DL-SCH). The transmitter model includes PDSCH demodulation reference signals (DM-RS) and PDSCH phase tracking reference signals (PT-RS). The example supports NTN narrowband and NTN tapped delay line (TDL) propagation channels.

Introduction

This example measures the PDSCH throughput of a 5G link, as defined by the 3GPP NR standard [1], [2], [3], [4].

The example models these 5G NR features:

  • DL-SCH transport channel coding

  • Multiple codewords, dependent on the number of layers

  • PDSCH, PDSCH DM-RS, and PDSCH PT-RS generation

  • Variable subcarrier spacing and frame numerologies

  • Normal and extended cyclic prefix

  • NTN narrowband and NTN TDL propagation channel models

Other features of the simulation are:

  • PDSCH precoding using singular value decomposition (SVD)

  • Cyclic prefix orthogonal frequency division multiplexing (CP-OFDM) modulation

  • Slot-wise and non-slot-wise PDSCH and DM-RS mapping

  • Timing synchronization and channel estimation

  • A single bandwidth part (BWP) across the whole carrier

The figure shows the implemented processing chain. For clarity, DM-RS and PT-RS generation are omitted.

modeling_NTN_Channel_image_2.jpg

For a more detailed explanation of the steps implemented in this example, see Model 5G NR Communication Links (5G Toolbox) and DL-SCH and PDSCH Transmit and Receive Processing Chain (5G Toolbox).

This example supports wideband and subband precoding. You determine the precoding matrix by using SVD and averaging the channel estimate across all PDSCH PRBs either in the allocation (for wideband precoding) or in the subband.

To reduce the total simulation time, you can use Parallel Computing Toolbox™ to execute the SNR points of the SNR loop in parallel.

Simulation Length and SNR Points

Set the length of the simulation in terms of the number of 10 ms frames. By default, the example uses 2 frames, but, a large number of 10 ms frames is necessary to produce meaningful throughput results. Set the SNR points to simulate. The SNR for each layer is defined per resource element (RE), and it includes the effect of signal and noise across all antennas. For an explanation of the SNR definition that this example uses, see SNR Definition Used in Link Simulations (5G Toolbox).

simParameters = struct;         % Clear simParameters variable to contain 
                                % all key simulation parameters
simParameters.NFrames = 2;      % Number of 10 ms frames
simParameters.SNRIn = 5:2:15;   % SNR range (dB)

Set the displaySimulationInformation variable to true to display information about the throughput simulation at each SNR point.

displaySimulationInformation = true;

Doppler Compensation Configuration

The example supports two Doppler compensation configurations: one at transmitter end and the other at receiver end. For compensation at transmitter end, use the DopplerPreCompensator field. Set the DopplerPreCompensator field to true to account for Doppler due to satellite movement by applying Doppler pre-compensation to the transmitted waveform. For compensation at receiver end, use the RxDopplerCompensator field. Set the RxDopplerCompensator field to true to estimate and compensate the Doppler shift for the received waveform.

simParameters.DopplerPreCompensator = true;
simParameters.RxDopplerCompensator = false;

Carrier and PDSCH Configuration

Set the key parameters of the simulation. These parameters include:

  • Bandwidth in resource blocks (RBs)

  • Subcarrier spacing (SCS) in kHz: 15, 30, 60, or 120

  • Cyclic prefix length (CP): normal or extended

  • Cell identity

  • Number of transmit and receive antennas

Also create a substructure containing the DL-SCH and PDSCH parameters, including:

  • Target code rate

  • Allocated resource blocks (PRBSet)

  • Modulation scheme: 'QPSK', '16QAM', '64QAM', or '256QAM'

  • Number of layers

  • PDSCH mapping type

  • DM-RS configuration parameters

  • PT-RS configuration parameters

% Set waveform type and PDSCH numerology (SCS and CP type)
simParameters.Carrier = nrCarrierConfig;
simParameters.Carrier.SubcarrierSpacing = 30;
simParameters.Carrier.CyclicPrefix = 'Normal';
% Bandwidth in number of RBs (11 RBs at 30 kHz SCS for 5 MHz bandwidth)
simParameters.Carrier.NSizeGrid = 11;
% Physical layer cell identity
simParameters.Carrier.NCellID = 1;

% PDSCH/DL-SCH parameters
% This PDSCH definition is the basis for all PDSCH transmissions in the
% throughput simulation
simParameters.PDSCH = nrPDSCHConfig;
% This structure is to hold additional simulation parameters for the DL-SCH
% and PDSCH
simParameters.PDSCHExtension = struct();

% Define PDSCH time-frequency resource allocation per slot to be full grid
% (single full grid BWP)
% PDSCH PRB allocation
simParameters.PDSCH.PRBSet = 0:simParameters.Carrier.NSizeGrid-1;
% Starting symbol and number of symbols of each PDSCH allocation
simParameters.PDSCH.SymbolAllocation = [0,simParameters.Carrier.SymbolsPerSlot];
simParameters.PDSCH.MappingType = 'A';

% Scrambling identifiers
simParameters.PDSCH.NID = simParameters.Carrier.NCellID;
simParameters.PDSCH.RNTI = 1;

% PDSCH resource block mapping (TS 38.211 Section 7.3.1.6)
simParameters.PDSCH.VRBToPRBInterleaving = 0;
simParameters.PDSCH.VRBBundleSize = 4;

% Define the number of transmission layers to be used
simParameters.PDSCH.NumLayers = 1;

% Define codeword modulation and target coding rate
% The number of codewords is directly dependent on the number of layers so
% ensure that layers are set first before getting the codeword number
if simParameters.PDSCH.NumCodewords > 1
    % Multicodeword transmission (when number of layers being > 4)
    simParameters.PDSCH.Modulation = {'16QAM','16QAM'};
    % Code rate used to calculate transport block sizes
    simParameters.PDSCHExtension.TargetCodeRate = [490 490]/1024;
else
    simParameters.PDSCH.Modulation = '16QAM';
    % Code rate used to calculate transport block size
    simParameters.PDSCHExtension.TargetCodeRate = 490/1024;
end

% DM-RS and antenna port configuration (TS 38.211 Section 7.4.1.1)
simParameters.PDSCH.DMRS.DMRSPortSet = 0:simParameters.PDSCH.NumLayers-1;
simParameters.PDSCH.DMRS.DMRSTypeAPosition = 2;
simParameters.PDSCH.DMRS.DMRSLength = 1;
simParameters.PDSCH.DMRS.DMRSAdditionalPosition = 2;
simParameters.PDSCH.DMRS.DMRSConfigurationType = 2;
simParameters.PDSCH.DMRS.NumCDMGroupsWithoutData = 1;
simParameters.PDSCH.DMRS.NIDNSCID = 1;
simParameters.PDSCH.DMRS.NSCID = 0;

% PT-RS configuration (TS 38.211 Section 7.4.1.2)
simParameters.PDSCH.EnablePTRS = 0;
simParameters.PDSCH.PTRS.TimeDensity = 1;
simParameters.PDSCH.PTRS.FrequencyDensity = 2;
simParameters.PDSCH.PTRS.REOffset = '00';
% PT-RS antenna port, subset of DM-RS port set. Empty corresponds to lowest
% DM-RS port number
simParameters.PDSCH.PTRS.PTRSPortSet = [];

% Reserved PRB patterns, if required (for CORESETs, forward compatibility etc)
simParameters.PDSCH.ReservedPRB{1}.SymbolSet = [];   % Reserved PDSCH symbols
simParameters.PDSCH.ReservedPRB{1}.PRBSet = [];      % Reserved PDSCH PRBs
simParameters.PDSCH.ReservedPRB{1}.Period = [];      % Periodicity of reserved resources

% Additional simulation and DL-SCH related parameters
% PDSCH PRB bundling (TS 38.214 Section 5.1.2.3)
simParameters.PDSCHExtension.PRGBundleSize = [];     % 2, 4, or [] to signify "wideband"
% Rate matching or transport block size (TBS) parameters
% Set PDSCH rate matching overhead for TBS (Xoh) to 6 when PT-RS is enabled, otherwise 0
simParameters.PDSCHExtension.XOverhead = 6*simParameters.PDSCH.EnablePTRS;
% Hybrid automatic repeat request (HARQ) parameters
% Number of parallel HARQ processes to use
simParameters.PDSCHExtension.NHARQProcesses = 1;
% Enable retransmissions for each process, using redundancy version (RV) sequence [0,2,3,1]
simParameters.PDSCHExtension.EnableHARQ = false;
% LDPC decoder parameters
% Available algorithms: 'Belief propagation', 'Layered belief propagation',
%                       'Normalized min-sum', 'Offset min-sum'
simParameters.PDSCHExtension.LDPCDecodingAlgorithm = 'Normalized min-sum';
simParameters.PDSCHExtension.MaximumLDPCIterationCount = 6;

% Define the overall transmission antenna geometry at end-points
% For NTN narrowband channel, only single-input-single-output (SISO)
% transmission is allowed
% Number of PDSCH transmission antennas (1,2,4,8,16,32,64,128,256,512,1024) >= NumLayers
simParameters.NumTransmitAntennas = 1;
if simParameters.PDSCH.NumCodewords > 1 % Multi-codeword transmission
    % Number of UE receive antennas (even number >= NumLayers)
    simParameters.NumReceiveAntennas = 8;
else
    % Number of UE receive antennas (1 or even number >= NumLayers)
    simParameters.NumReceiveAntennas = 1;
end

Get information about the baseband waveform after the OFDM modulation step.

waveformInfo = nrOFDMInfo(simParameters.Carrier);

Propagation Channel Model Construction

Create the channel model object for the simulation. Both the NTN narrowband and NTN TDL channel models are supported [5], [6]. For more information on how to model NTN narrowband and NTN TDL channels, see Model NR NTN Channel.

% Define the general NTN propagation channel parameters
% Set the NTN channel type to Narrowband for an NTN narrowband channel and
% set the NTN channel type to TDL for a NTN TDL channel.
simParameters.NTNChannelType = 'Narrowband';

% Set the parameters common to both NTN narrowband and NTN TDL channels
simParameters.CarrierFrequency = 2e9;               % Carrier frequency (in Hz)
simParameters.ElevationAngle = 50;                  % Elevation angle (in degrees)
simParameters.MobileSpeed = 3*1000/3600;            % Speed of mobile terminal (in m/s)
simParameters.SatelliteAltitude = 600000;           % Satellite altitude (in m)
simParameters.SatelliteSpeed = 7562.2;              % Satellite speed (in m/s)
simParameters.SampleRate = waveformInfo.SampleRate;
simParameters.RandomStream = 'mt19937ar with seed';
simParameters.Seed = 73;

% Set the following fields for NTN narrowband channel
if strcmpi(simParameters.NTNChannelType,'Narrowband')
    simParameters.Environment = 'Urban';
    simParameters.AzimuthOrientation = 0;
end

% Set the following fields for NTN TDL channel
if strcmpi(simParameters.NTNChannelType,'TDL')
    simParameters.DelayProfile = 'NTN-TDL-A';
    simParameters.DelaySpread = 30e-9;
end

% Cross-check the PDSCH layering against the channel geometry 
validateNumLayers(simParameters);

% Set up the NTN channel
ntnChan = HelperSetupNTNChannel(simParameters);

% Get the maximum number of delayed samples due to a channel multipath
% component. The maximum number of delayed samples is calculated from the
% channel path with the maximum delay and the implementation delay of the
% channel filter. This number of delay samples is required later to flush
% the channel filter to obtain the received signal.
chInfo = info(ntnChan.BaseChannel);
maxChDelay = ceil(max(chInfo.PathDelays*ntnChan.BaseChannel.SampleRate)) + ...
    chInfo.ChannelFilterDelay;

Processing Loop

To determine the throughput at each SNR point, analyze the PDSCH data for each transmission instance using these steps.

  1. Generate the transport block — Get the transport block size for each codeword depending on the PDSCH configuration. Generate new transport block(s) for each transmission depending on the transmission status for given HARQ process.

  2. Generate the resource grid — The nrDLSCH (5G Toolbox) System object performs transport channel coding. The object operates on the input transport block. The nrPDSCH (5G Toolbox) function modulates the encoded data bits. Apply an implementation-specific multiple-input-multiple-output (MIMO) precoding to the modulated symbols. Map these modulated symbols along with reference signal to the resource grid.

  3. Generate the waveform — The nrOFDMModulate (5G Toolbox) function OFDM-modulates the generated grid to get the time-domain waveform.

  4. Apply Doppler pre-compensation — Apply the Doppler shift due to satellite movement to the generated waveform to pre-compensate the channel induced satellite Doppler shift.

  5. Model and apply a noisy channel — Pass the generated waveform through an NTN narrowband or NTN TDL fading channel to get the faded waveform. Add additive white Gaussian noise (AWGN) to the faded waveform. Define the SNR for each layer based on the RE and receive antenna. For an explanation of the SNR definition that this example uses, see SNR Definition Used in Link Simulations (5G Toolbox).

  6. Apply Doppler compensation — Estimate the Doppler shift in the received waveform and compensate the Doppler shift.

  7. Perform synchronization and OFDM demodulation — For timing synchronization, the received waveform is correlated with the PDSCH DM-RS. The nrOFDMDemodulate (5G Toolbox) function then OFDM-demodulates the synchronized signal.

  8. Perform channel estimation — For channel estimation, PDSCH DM-RS is used.

  9. Perform equalization and CPE compensation — The nrEqualizeMMSE (5G Toolbox) function equalizes the received PDSCH REs. Use the PT-RS symbols to estimate the common phase error (CPE) and then correct the error in each OFDM symbol within the range of the reference PT-RS OFDM symbols.

  10. Calculate precoding matrix — Use SVD to generate the precoding matrix W for the next transmission.

  11. Decode the PDSCH — Demodulate and descramble the equalized PDSCH symbols, along with a noise estimate using the nrPDSCHDecode (5G Toolbox) function to obtain an estimate of the received codewords.

  12. Decode the DL-SCH — Pass the decoded soft bits through the nrDLSCHDecoder (5G Toolbox) System object. The object decodes the codeword and returns the block cyclic redundancy check (CRC) error. Update the HARQ process with the CRC error. This example determines the throughput of the PDSCH link using the CRC error.

% Array to store the maximum throughput for all SNR points
maxThroughput = zeros(length(simParameters.SNRIn),1); 
% Array to store the simulation throughput for all SNR points
simThroughput = zeros(length(simParameters.SNRIn),1);

% Set up RV sequence for all HARQ processes
if simParameters.PDSCHExtension.EnableHARQ
    % In the final report of RAN WG1 meeting #91 (R1-1719301), it was
    % observed in R1-1717405 that if performance is the priority, [0 2 3 1]
    % should be used. If self-decodability is the priority, it should be
    % taken into account that the upper limit of the code rate at which
    % each RV is self-decodable is in the following order: 0>3>2>1
    rvSeq = [0 2 3 1];
else
    % In case of HARQ disabled, RV is set to 0
    rvSeq = 0;
end

% Create DL-SCH encoder System object to perform transport channel encoding
encodeDLSCH = nrDLSCH;
encodeDLSCH.MultipleHARQProcesses = true;
encodeDLSCH.TargetCodeRate = simParameters.PDSCHExtension.TargetCodeRate;

% Create DL-SCH decoder System object to perform transport channel decoding
decodeDLSCH = nrDLSCHDecoder;
decodeDLSCH.MultipleHARQProcesses = true;
decodeDLSCH.TargetCodeRate = simParameters.PDSCHExtension.TargetCodeRate;
decodeDLSCH.LDPCDecodingAlgorithm = simParameters.PDSCHExtension.LDPCDecodingAlgorithm;
decodeDLSCH.MaximumLDPCIterationCount = ...
    simParameters.PDSCHExtension.MaximumLDPCIterationCount;

% Set a threshold value to detect the valid OFDM symbol boundary. For a
% SISO case, a threshold of 0.48 can be used to have probability of
% incorrect boundary detection around 0.01. Use 0 to avoid thresholding
% logic.
dtxThresold = 0.48;

% Use an offset to account for the common delay. The example, by default,
% does not introduce any common delay and only passes through the channel.
sampleDelayOffset = 0; % Number of samples

% Set this flag to true, to use the shift value estimated in first slot
% directly for the consecutive slots. When set to false, the shift is
% calculated for each slot, considering the range of shift values to be
% whole cyclic prefix length. This is used in the estimation of integer
% Doppler shift.
usePreviousShift = false;

% Processing loop
for snrIdx = 1:numel(simParameters.SNRIn)      % Comment out for parallel computing
% parfor snrIdx = 1:numel(simParameters.SNRIn) % Uncomment for parallel computing
    % To reduce the total simulation time, you can execute this loop in
    % parallel by using Parallel Computing Toolbox features. Comment
    % out the for-loop statement and uncomment the parfor-loop statement.
    % If Parallel Computing Toolbox is not installed, parfor-loop defaults
    % to a for-loop statement. Because the parfor-loop iterations are
    % executed in parallel in a nondeterministic order, the simulation
    % information displayed for each SNR point can be intertwined. To
    % switch off the simulation information display, set the
    % displaySimulationInformation variable (defined earlier in this
    % example) to false.

    % Reset the random number generator so that each SNR point experiences
    % the same noise realization
    rng('default');

    % Make copies of the simulation-level parameter structures so that they
    % are not Parallel Computing Toolbox broadcast variables when using parfor
    simLocal = simParameters;
    waveinfoLocal = waveformInfo;
    ntnChanLocal = ntnChan;

    % Make copies of channel-level parameters to simplify subsequent
    % parameter referencing
    carrier = simLocal.Carrier;
    pdsch = simLocal.PDSCH;
    pdschextra = simLocal.PDSCHExtension;
    % Copy of the decoder handle to help Parallel Computing Toolbox
    % classification
    decodeDLSCHLocal = decodeDLSCH;
    decodeDLSCHLocal.reset();       % Reset decoder at the start of each SNR point
    pathFilters = [];
    thres = dtxThresold;
    sampleOffset = sampleDelayOffset;
    usePrevShift = usePreviousShift;

    % SNR value
    SNRdB = simLocal.SNRIn(snrIdx);

    % Specify the order in which we cycle through the HARQ process
    % identifiers
    harqSequence = 0:pdschextra.NHARQProcesses-1;

    % Initialize the state of all HARQ processes
    harqEntity = HARQEntity(harqSequence,rvSeq,pdsch.NumCodewords);

    % Reset the channel so that each SNR point experiences the same channel
    % realization
    reset(ntnChanLocal.BaseChannel);
    reset(ntnChanLocal.ChannelFilter);

    % Total number of slots in the simulation period
    NSlots = simLocal.NFrames*carrier.SlotsPerFrame;

    % Obtain a precoding matrix (wtx) to use in the transmission of the
    % first transport block
    [estChannelGrid,sampleTimes] = getInitialChannelEstimate(...
        carrier,simLocal.NumTransmitAntennas,ntnChanLocal);    
    newWtx = getPrecodingMatrix(carrier,pdsch,estChannelGrid,pdschextra.PRGBundleSize);

    % Timing offset, updated in every slot when the correlation is strong
    offset = 0;
    tDoppler = sampleTimes(end);
    shiftOut = 0;

    % Loop over the entire waveform length
    for nslot = 0:NSlots-1

        % Update carrier slot number to account for new slot transmission
        carrier.NSlot = nslot;

        % Calculate the transport block sizes for the transmission in the slot
        [pdschIndices,pdschIndicesInfo] = nrPDSCHIndices(carrier,pdsch);
        trBlkSizes = nrTBS(pdsch.Modulation,pdsch.NumLayers,numel(pdsch.PRBSet),...
            pdschIndicesInfo.NREPerPRB,pdschextra.TargetCodeRate,pdschextra.XOverhead);

        % Set transport block depending on the HARQ process
        for cwIdx = 1:pdsch.NumCodewords
            % Create a new DL-SCH transport block for new data in the
            % current process
            if harqEntity.NewData(cwIdx) 
                trBlk = randi([0 1],trBlkSizes(cwIdx),1);
                setTransportBlock(encodeDLSCH,trBlk,cwIdx-1,harqEntity.HARQProcessID);
                % Flush decoder soft buffer explicitly for any new data
                % because of previous RV sequence time out
                if harqEntity.SequenceTimeout(cwIdx)
                    resetSoftBuffer(decodeDLSCHLocal,cwIdx-1,harqEntity.HARQProcessID);
                end
            end
        end

        % Encode the DL-SCH transport blocks
        codedTrBlocks = encodeDLSCH(pdsch.Modulation,pdsch.NumLayers, ...
            pdschIndicesInfo.G,harqEntity.RedundancyVersion,harqEntity.HARQProcessID);

        % Get precoding matrix (wtx) calculated in previous slot
        wtx = newWtx;

        % Perform PDSCH modulation
        pdschSymbols = nrPDSCH(carrier,pdsch,codedTrBlocks);

        % Create resource grid associated with PDSCH transmission antennas
        pdschGrid = nrResourceGrid(carrier,simLocal.NumTransmitAntennas);

        % Perform implementation-specific PDSCH MIMO precoding and mapping
        [pdschAntSymbols,pdschAntIndices] = hPRGPrecode(size(pdschGrid), ...
            carrier.NStartGrid,pdschSymbols,pdschIndices,wtx);
        pdschGrid(pdschAntIndices) = pdschAntSymbols;

        % Perform implementation-specific PDSCH DM-RS MIMO precoding and
        % mapping
        dmrsSymbols = nrPDSCHDMRS(carrier,pdsch);
        dmrsIndices = nrPDSCHDMRSIndices(carrier,pdsch);
        [dmrsAntSymbols,dmrsAntIndices] = hPRGPrecode(size(pdschGrid), ...
            carrier.NStartGrid,dmrsSymbols,dmrsIndices,wtx);        
        pdschGrid(dmrsAntIndices) = dmrsAntSymbols;

        % Perform implementation-specific PDSCH PT-RS MIMO precoding and
        % mapping
        ptrsSymbols = nrPDSCHPTRS(carrier,pdsch);
        ptrsIndices = nrPDSCHPTRSIndices(carrier,pdsch);
        [ptrsAntSymbols,ptrsAntIndices] = hPRGPrecode(size(pdschGrid), ...
            carrier.NStartGrid,ptrsSymbols,ptrsIndices,wtx);
        pdschGrid(ptrsAntIndices) = ptrsAntSymbols;        

        % Perform OFDM modulation
        txWaveform = nrOFDMModulate(carrier,pdschGrid);

        % Pass data through the channel model. Append zeros at the end of
        % the transmitted waveform to flush the channel content. These
        % zeros take into account any delay introduced in the channel. This
        % delay is a combination of the multipath delay and implementation
        % delay. This value can change depending on the sampling rate,
        % delay profile, and delay spread. Also apply Doppler
        % pre-compensation.
        txWaveform = [txWaveform; zeros(maxChDelay,size(txWaveform,2))]; %#ok<AGROW>
        [txWaveform,tDoppler] = compensateDopplerShift(...
            txWaveform,ntnChanLocal.BaseChannel.SampleRate, ...
            ntnChanLocal.SatelliteDopplerShift, ...
            simLocal.DopplerPreCompensator,tDoppler);
        [rxWaveform,pathGains] = HelperGenerateNTNChannel(...
            ntnChanLocal,txWaveform);

        % Add AWGN to the received time-domain waveform. Normalize the
        % noise power by the size of the inverse fast Fourier transform
        % (IFFT) used in OFDM modulation, because the OFDM modulator
        % applies this normalization to the transmitted waveform. Also,
        % normalize the noise power by the number of receive antennas,
        % because the default behavior of the channel model is to apply
        % this normalization to the received waveform.
        SNR = 10^(SNRdB/10);
        noiseScaling = ...
            1/sqrt(2.0*simLocal.NumReceiveAntennas*double(waveinfoLocal.Nfft)*SNR);
        noise = noiseScaling*complex(randn(size(rxWaveform)),randn(size(rxWaveform)));
        rxWaveform = rxWaveform + noise;

        % Perform fractional Doppler frequency shift estimation and
        % compensation. Use the cyclic prefix in the OFDM waveform to
        % compute the fractional Doppler shift.
        [fractionalDopplerShift,detFlag] = estimateFractionalDopplerShift( ...
            rxWaveform,carrier.SubcarrierSpacing,waveinfoLocal.Nfft, ...
            waveinfoLocal.CyclicPrefixLengths(2),thres, ...
            simLocal.RxDopplerCompensator);
        rxWaveform = compensateDopplerShift(rxWaveform,waveinfoLocal.SampleRate, ...
            fractionalDopplerShift,simLocal.RxDopplerCompensator);

        % Perform integer Doppler frequency shift estimation and
        % compensation. Use the demodulation reference signals to compute
        % the integer Doppler shift.
        [integerDopplerShift,shiftOut] = estimateIntegerDopplerShift(carrier,rxWaveform, ...
            dmrsIndices,dmrsSymbols,sampleOffset,usePrevShift,shiftOut-sampleOffset, ...
            (simLocal.RxDopplerCompensator && detFlag));
        rxWaveform = compensateDopplerShift(rxWaveform,waveinfoLocal.SampleRate, ...
            integerDopplerShift,simLocal.RxDopplerCompensator);

        % For timing synchronization, correlate the received waveform with
        % the PDSCH DM-RS to give timing offset estimate t and correlation
        % magnitude mag. The function hSkipWeakTimingOffset is used to
        % update the receiver timing offset. If the correlation peak in mag
        % is weak, the current timing estimate t is ignored and the
        % previous estimate offset is used.
        [t,mag] = nrTimingEstimate(carrier,rxWaveform,dmrsIndices,dmrsSymbols); 
        offset = hSkipWeakTimingOffset(offset,t,mag);
        % Display a warning if the estimated timing offset exceeds the
        % maximum channel delay
        if offset > maxChDelay
            warning(['Estimated timing offset (%d) is greater than the maximum' ...
                ' channel delay (%d). This might result in a decoding failure.' ...
                ' The estimated timing offset might be greater than the maximum' ...
                ' channel delay because of low SNR or not enough DM-RS symbols to' ...
                ' synchronize successfully.'],offset,maxChDelay);
        end
        rxWaveform = rxWaveform(1+offset:end,:);

        % Perform OFDM demodulation on the received data to recreate the
        % resource grid. Include zero padding in the event that practical
        % synchronization results in an incomplete slot being demodulated.
        rxGrid = nrOFDMDemodulate(carrier,rxWaveform);
        [K,L,R] = size(rxGrid);
        if (L < carrier.SymbolsPerSlot)
            rxGrid = cat(2,rxGrid,zeros(K,carrier.SymbolsPerSlot-L,R));
        end

        % Perform least squares channel estimation between the received
        % grid and each transmission layer, using the PDSCH DM-RS for each
        % layer. This channel estimate includes the effect of transmitter
        % precoding.
        [estChannelGrid,noiseEst] = nrChannelEstimate(carrier,rxGrid,...
            dmrsIndices,dmrsSymbols,'CDMLengths',pdsch.DMRS.CDMLengths); 

        % Get PDSCH REs from the received grid and estimated channel grid
        [pdschRx,pdschHest] = nrExtractResources(...
            pdschIndices,rxGrid,estChannelGrid);

        % Remove precoding from estChannelGrid prior to precoding
        % matrix calculation
        estChannelGridPorts = precodeChannelEstimate(...
            carrier,estChannelGrid,conj(wtx));

        % Get precoding matrix for next slot
        newWtx = getPrecodingMatrix(...
            carrier,pdsch,estChannelGridPorts,pdschextra.PRGBundleSize);

        % Perform equalization
        [pdschEq,csi] = nrEqualizeMMSE(pdschRx,pdschHest,noiseEst);

        % Common phase error (CPE) compensation
        if ~isempty(ptrsIndices)
            % Initialize temporary grid to store equalized symbols
            tempGrid = nrResourceGrid(carrier,pdsch.NumLayers);

            % Extract PT-RS symbols from received grid and estimated
            % channel grid
            [ptrsRx,ptrsHest,~,~,~,ptrsLayerIndices] = ...
                nrExtractResources(ptrsIndices,rxGrid,estChannelGrid,tempGrid);

            % Equalize PT-RS symbols and map them to tempGrid
            ptrsEq = nrEqualizeMMSE(ptrsRx,ptrsHest,noiseEst);
            tempGrid(ptrsLayerIndices) = ptrsEq;

            % Estimate the residual channel at the PT-RS locations in
            % tempGrid
            cpe = nrChannelEstimate(tempGrid,ptrsIndices,ptrsSymbols);

            % Sum estimates across subcarriers, receive antennas, and
            % layers. Then, get the CPE by taking the angle of the
            % resultant sum
            cpe = angle(sum(cpe,[1 3 4]));

            % Map the equalized PDSCH symbols to tempGrid
            tempGrid(pdschIndices) = pdschEq;

            % Correct CPE in each OFDM symbol within the range of reference
            % PT-RS OFDM symbols
            symLoc = ...
                pdschIndicesInfo.PTRSSymbolSet(1)+1:pdschIndicesInfo.PTRSSymbolSet(end)+1;
            tempGrid(:,symLoc,:) = tempGrid(:,symLoc,:).*exp(-1i*cpe(symLoc));

            % Extract PDSCH symbols
            pdschEq = tempGrid(pdschIndices);
        end

        % Decode PDSCH symbols
        [dlschLLRs,rxSymbols] = nrPDSCHDecode(carrier,pdsch,pdschEq,noiseEst);

        % Scale the decoded log-likelihood ratios (LLRs) by channel state
        % information (CSI)
        csi = nrLayerDemap(csi); % CSI layer demapping
        for cwIdx = 1:pdsch.NumCodewords
            Qm = length(dlschLLRs{cwIdx})/length(rxSymbols{cwIdx}); % Bits per symbol
            csi{cwIdx} = repmat(csi{cwIdx}.',Qm,1);                 % Expand by each bit
                                                                    % per symbol
            dlschLLRs{cwIdx} = dlschLLRs{cwIdx} .* csi{cwIdx}(:);   % Scale by CSI
        end

        % Decode the DL-SCH transport channel
        decodeDLSCHLocal.TransportBlockLength = trBlkSizes;
        [decbits,blkerr] = decodeDLSCHLocal(dlschLLRs,pdsch.Modulation,...
            pdsch.NumLayers,harqEntity.RedundancyVersion,harqEntity.HARQProcessID);

        % Store values to calculate throughput
        simThroughput(snrIdx) = simThroughput(snrIdx) + sum(~blkerr .* trBlkSizes);
        maxThroughput(snrIdx) = maxThroughput(snrIdx) + sum(trBlkSizes);

        % Update current process with CRC error and advance to next process
        updateAndAdvance(harqEntity,blkerr,trBlkSizes,pdschIndicesInfo.G);
     end

    % Display the results
    if displaySimulationInformation == 1
        fprintf('\nThroughput(Mbps) for %d frame(s) at SNR %d dB: %.4f\n',...
            simLocal.NFrames,SNRdB,1e-6*simThroughput(snrIdx)/(simLocal.NFrames*10e-3));
        fprintf('Throughput(%%) for %d frame(s) at SNR %d dB: %.4f\n',...
            simLocal.NFrames,SNRdB,simThroughput(snrIdx)*100/maxThroughput(snrIdx));
    end

end
Throughput(Mbps) for 2 frame(s) at SNR 5 dB: 0.0000
Throughput(%) for 2 frame(s) at SNR 5 dB: 0.0000
Throughput(Mbps) for 2 frame(s) at SNR 7 dB: 0.0000
Throughput(%) for 2 frame(s) at SNR 7 dB: 0.0000
Throughput(Mbps) for 2 frame(s) at SNR 9 dB: 6.2308
Throughput(%) for 2 frame(s) at SNR 9 dB: 92.5000
Throughput(Mbps) for 2 frame(s) at SNR 11 dB: 6.7360
Throughput(%) for 2 frame(s) at SNR 11 dB: 100.0000
Throughput(Mbps) for 2 frame(s) at SNR 13 dB: 6.7360
Throughput(%) for 2 frame(s) at SNR 13 dB: 100.0000
Throughput(Mbps) for 2 frame(s) at SNR 15 dB: 6.7360
Throughput(%) for 2 frame(s) at SNR 15 dB: 100.0000

Results

Display the measured throughput, which is the percentage of the maximum possible throughput of the link given the available resources for data transmission.

figure;
plot(simParameters.SNRIn,simThroughput*100./maxThroughput,'o-.')
xlabel('SNR (dB)'); ylabel('Throughput (%)'); grid on;
title(sprintf('NTN %s (%dx%d) / NRB=%d / SCS=%dkHz', ...
              simParameters.NTNChannelType,simParameters.NumTransmitAntennas, ...
              simParameters.NumReceiveAntennas, simParameters.Carrier.NSizeGrid,...
              simParameters.Carrier.SubcarrierSpacing));

Figure contains an axes object. The axes object with title NTN Narrowband (1x1) / NRB=11 / SCS=30kHz contains an object of type line.

% Bundle key parameters and results into a combined structure for recording
simResults.simParameters = simParameters;
simResults.simThroughput = simThroughput;
simResults.maxThroughput = maxThroughput;

This next figure shows the throughput results obtained by simulating 1000 frames (NFrames = 1000, SNRIn = 5:2:15) for a carrier with a 30 kHz SCS occupying a 5 MHz transmission bandwidth. The simulation setup includes the default carrier and PDSCH configuration with an NTN narrowband channel. The line corresponding to the Doppler Pre-compensator is achieved by setting the DopplerPreCompensator field to true and RxDopplerCompensator field to false. The line corresponding to the Rx Doppler Compensator is achieved by setting the DopplerPreCompensator field to false and RxDopplerCompensator field to true.

example_Throughput_Figure_U_2.jpg

Further Exploration

You can use this example to further explore these options:

  • To analyze the throughput at each SNR point for a different satellite orbit, vary the satellite altitude and satellite speed.

  • To observe the link performance without any Doppler compensation techniques, set DopplerPreCompensator and RxDopplerCompensator fields to false.

  • To observe the link performance in case of no Doppler pre-compensation and using receiver techniques to compensate Doppler shift, set the DopplerPreCompensator field to false and RxDopplerCompensator field to true.

  • To check the throughput performance of different scenarios, change the carrier numerology and the number of transmit and receive antennas, and set the channel model type to TDL.

  • To compare the throughput performance in an NTN and terrestrial network, use the nrTDLChannel (5G Toolbox) and the nrCDLChannel (5G Toolbox) channel objects as shown in NR PDSCH Throughput (5G Toolbox).

Appendix

The example uses these helper functions:

Selected Bibliography

[1] 3GPP TS 38.211. "NR; Physical channels and modulation." 3rd Generation Partnership Project; Technical Specification Group Radio Access Network.

[2] 3GPP TS 38.212. "NR; Multiplexing and channel coding." 3rd Generation Partnership Project; Technical Specification Group Radio Access Network.

[3] 3GPP TS 38.213. "NR; Physical layer procedures for control." 3rd Generation Partnership Project; Technical Specification Group Radio Access Network.

[4] 3GPP TS 38.214. "NR; Physical layer procedures for data." 3rd Generation Partnership Project; Technical Specification Group Radio Access Network.

[5] 3GPP TR 38.901. "Study on channel model for frequencies from 0.5 to 100 GHz." 3rd Generation Partnership Project; Technical Specification Group Radio Access Network.

[6] 3GPP TR 38.811. "Study on new radio (NR) to support non-terrestrial networks." 3rd Generation Partnership Project; Technical Specification Group Radio Access Network.

[7] 3GPP TR 38.821. "Solutions for NR to support non-terrestrial networks (NTN)." 3rd Generation Partnership Project; Technical Specification Group Radio Access Network.

[8] ITU-R Recommendation P.681-11 (08/2019). "Propagation data required for the design systems in the land mobile-satellite service." P Series; Radio wave propagation.

Local Functions

function validateNumLayers(simParameters)
% Validate the number of layers, relative to the antenna geometry

    numlayers = simParameters.PDSCH.NumLayers;
    ntxants = simParameters.NumTransmitAntennas;
    nrxants = simParameters.NumReceiveAntennas;

    if contains(simParameters.NTNChannelType,'Narrowband','IgnoreCase',true)
        if (ntxants ~= 1) || (nrxants ~= 1)
            error(['For NTN narrowband channel, ' ... 
                'the number of transmit and receive antennas must be 1.']);
        end
    end

    antennaDescription = sprintf(...
        'min(NumTransmitAntennas,NumReceiveAntennas) = min(%d,%d) = %d', ...
        ntxants,nrxants,min(ntxants,nrxants));
    if numlayers > min(ntxants,nrxants)
        error('The number of layers (%d) must satisfy NumLayers <= %s', ...
            numlayers,antennaDescription);
    end

    % Display a warning if the maximum possible rank of the channel equals
    % the number of layers
    if (numlayers > 2) && (numlayers == min(ntxants,nrxants))
        warning(['The maximum possible rank of the channel, given by %s, is equal to' ...
            ' NumLayers (%d). This may result in a decoding failure under some channel' ...
            ' conditions. Try decreasing the number of layers or increasing the channel' ...
            ' rank (use more transmit or receive antennas).'],antennaDescription, ...
            numlayers); %#ok<SPWRN>
    end

end

function [estChannelGrid,sampleTimes] = getInitialChannelEstimate(...
    carrier,nTxAnts,ntnChan)
% Obtain channel estimate before first transmission. Use this function to
% obtain a precoding matrix for the first slot.

    ofdmInfo = nrOFDMInfo(carrier);

    chInfo = info(ntnChan.BaseChannel);
    maxChDelay = ceil(max(chInfo.PathDelays*ntnChan.BaseChannel.SampleRate)) ...
        + chInfo.ChannelFilterDelay;

    % Temporary waveform (only needed for the sizes)
    tmpWaveform = zeros(...
        (ofdmInfo.SampleRate/1000/carrier.SlotsPerSubframe)+maxChDelay,nTxAnts);

    % Filter through channel and get the path gains and path filters
    [~,pathGains,sampleTimes] = HelperGenerateNTNChannel(...
        ntnChan,tmpWaveform);
    chanInfo = info(ntnChan.ChannelFilter);
    pathFilters = chanInfo.ChannelFilterCoefficients.';

    % Perfect timing synchronization
    offset = nrPerfectTimingEstimate(pathGains,pathFilters);

    % Perfect channel estimate
    estChannelGrid = nrPerfectChannelEstimate(...
        carrier,pathGains,pathFilters,offset,sampleTimes);

end

function wtx = getPrecodingMatrix(carrier,pdsch,hestGrid,prgbundlesize)
% Calculate precoding matrices for all precoding resource block groups
% (PRGs) in the carrier that overlap with the PDSCH allocation

    % Maximum common resource block (CRB) addressed by carrier grid
    maxCRB = carrier.NStartGrid + carrier.NSizeGrid - 1;

    % PRG size
    if nargin==4 && ~isempty(prgbundlesize)
        Pd_BWP = prgbundlesize;
    else
        Pd_BWP = maxCRB + 1;
    end

    % PRG numbers (1-based) for each RB in the carrier grid
    NPRG = ceil((maxCRB + 1) / Pd_BWP);
    prgset = repmat((1:NPRG),Pd_BWP,1);
    prgset = prgset(carrier.NStartGrid + (1:carrier.NSizeGrid).');

    [~,~,R,P] = size(hestGrid);
    wtx = zeros([pdsch.NumLayers P NPRG]);
    for i = 1:NPRG

        % Subcarrier indices within current PRG and within the PDSCH
        % allocation
        thisPRG = find(prgset==i) - 1;
        thisPRG = intersect(thisPRG,pdsch.PRBSet(:) + carrier.NStartGrid,'rows');
        prgSc = (1:12)' + 12*thisPRG';
        prgSc = prgSc(:);

        if (~isempty(prgSc))

            % Average channel estimate in PRG
            estAllocGrid = hestGrid(prgSc,:,:,:);
            Hest = permute(mean(reshape(estAllocGrid,[],R,P)),[2 3 1]);

            % SVD decomposition
            [~,~,V] = svd(Hest);
            wtx(:,:,i) = V(:,1:pdsch.NumLayers).';

        end

    end

    wtx = wtx / sqrt(pdsch.NumLayers); % Normalize by NumLayers

end

function estChannelGrid = precodeChannelEstimate(carrier,estChannelGrid,W)
% Apply precoding matrix W to the last dimension of the channel estimate

    [K,L,R,P] = size(estChannelGrid);
    estChannelGrid = reshape(estChannelGrid,[K*L R P]);
    estChannelGrid = hPRGPrecode([K L R P],carrier.NStartGrid,estChannelGrid, ...
        reshape(1:numel(estChannelGrid),[K*L R P]),W);
    estChannelGrid = reshape(estChannelGrid,K,L,R,[]);

end

function [loc,wMovSum,pho,bestAnt] = detectOFDMSymbolBoundary(rxWave,nFFT,cpLen,thres)
% Detect OFDM symbol boundary by calculating correlation of cyclic prefix

    % Capture the dimensions of received waveform
    [NSamples,R] = size(rxWave);

    % Append zeros of length nFFT across each receive antenna to the
    % received waveform
    waveformZeroPadded = [rxWave;zeros(nFFT,R)];

    % Get the portion of waveform till the last nFFT samples
    wavePortion1 = waveformZeroPadded(1:end-nFFT,:);

    % Get the portion of waveform delayed by nFFT
    wavePortion2 = waveformZeroPadded(1+nFFT:end,:);

    % Get the energy of each sample in both the waveform portions
    eWavePortion1 = abs(wavePortion1).^2;
    eWavePortion2 = abs(wavePortion2).^2;

    % Initialize the temporary variables
    wMovSum = zeros([NSamples R]);
    wEnergyPortion1 = zeros([NSamples+cpLen-1 R]);
    wEnergyPortion2 = wEnergyPortion1;

    % Perform correlation for each sample with the sample delayed by nFFT
    waveCorr = wavePortion1.*conj(wavePortion2);
    % Calculate the moving sum value for each cpLen samples, across each
    % receive antenna
    oneVec = ones(cpLen,1);
    for i = 1:R
        wConv = conv(waveCorr(:,i),oneVec);
        wMovSum(:,i) = wConv(cpLen:end);
        wEnergyPortion1(:,i) = conv(eWavePortion1(:,i),oneVec);
        wEnergyPortion2(:,i) = conv(eWavePortion2(:,i),oneVec);
    end

    % Get the normalized correlation value for the moving sum matrix
    pho = abs(wMovSum)./ ...
        (eps+sqrt(wEnergyPortion1(cpLen:end,:).*wEnergyPortion2(cpLen:end,:)));

    % Detect the peak locations in each receive antenna based on the
    % threshold. These peak locations correspond to the starting location
    % of each OFDM symbol in the received waveform.
    loc = cell(R,1);
    m = zeros(R,1);
    phoFactor = ceil(NSamples/nFFT);
    phoExt = [pho; -1*ones(phoFactor*nFFT - NSamples,R)];
    for col = 1:R
        p1 = reshape(phoExt(:,i),[],phoFactor);
        [pks,locTemp] = max(p1);
        locTemp = locTemp + (0:phoFactor-1).*nFFT;
        indicesToConsider = pks>=thres;
        loc{col} = locTemp(indicesToConsider);
        m(col) = max(pks);
    end
    bestAnt = find(m == max(m));

end

function [out,detFlag] = estimateFractionalDopplerShift(rxWave,scs, ...
    nFFT,cpLen,thres,flag)
% Estimate the fractional Doppler shift using cyclic prefix

    if flag
        % Detect the OFDM boundary locations
        [loc,wMovSum,~,bestAnt] = detectOFDMSymbolBoundary(rxWave, ...
            nFFT,cpLen,thres);

        % Get the average correlation value at the peak locations for the
        % first receive antenna having maximum correlation value
        wSamples = nan(1,1);
        if ~isempty(loc{bestAnt(1)})
            wSamples(1) = mean(wMovSum(loc{bestAnt(1)},bestAnt(1)));
        end

        % Compute the fractional Doppler shift
        if ~all(isnan(wSamples))
            out = -(mean(angle(wSamples),'omitnan')*scs*1e3)/(2*pi);
            % Flag to indicate that there is at least one OFDM symbol
            % detected
            detFlag = 1;
        else
            out = 0;
            detFlag = 0;
        end
    else
        out = 0;
        detFlag = 0;
    end

end

function [out,shiftOut] = estimateIntegerDopplerShift(carrier,rx,refInd, ...
    refSym,sampleOffset,usePrevShift,shiftIn,flag)
% Estimate the integer Doppler shift using demodulation reference signal

    if flag
        % Get OFDM information
        ofdmInfo = nrOFDMInfo(carrier);
        K = carrier.NSizeGrid*12;                % Number of subcarriers     
        L = ofdmInfo.SymbolsPerSlot;             % Number of OFDM symbols in slot
        P = ceil(max(double(refInd(:))/(K*L)));  % Number of layers
        cpLen = ofdmInfo.CyclicPrefixLengths(1); % Highest cyclic prefix length

        % Range of shift values to be used in integer frequency offset
        % estimation
        shiftValues = sampleOffset + shiftIn;
        if ~(usePrevShift && (shiftIn > 0))
            % Update range of shift values such that whole cyclic prefix
            % length is covered
            shiftValues = sampleOffset + (1:cpLen);
        end

        % Initialize temporary variables
        shiftLen = length(shiftValues);
        maxValue = complex(zeros(shiftLen,1));
        binIndex = zeros(shiftLen,1);
        R = size(rx,2);
        rx1 = [rx; zeros(1*(mod(size(rx,1),2)),R)]; % Append zero, if required
        rxLen = size(rx1,1);
        x_wave = zeros([rxLen P]);

        % Generate reference waveform
        refGrid = complex(zeros([K L P]));
        refGrid(refInd) = refSym;
        refWave = nrOFDMModulate(carrier,refGrid,'Windowing',0);
        refWave = [refWave; zeros((rxLen-size(refWave,1)),P)];

        % Find the fast Fourier transform (FFT) bin corresponding to
        % maximum correlation value for each shift value
        for shiftIdx = 1:shiftLen
            % Use the waveform from the shift value and append zeros
            tmp = rx1(shiftValues(shiftIdx):end,:);
            rx = [tmp; zeros(rxLen-size(tmp,1),R)];

            % Compute the correlation of received waveform with reference
            % waveform across different layers and receive antennas
            for rIdx = 1:R
                for p = 1:P
                    x_wave(:,rIdx,p) = ...
                        rx(:,rIdx).*conj(refWave(1:length(rx(:,rIdx)),p));
                end
            end

            % Aggregate the correlated waveform across multiple ports and
            % compute energy of the resultant for each receive antenna
            x1 = sum(x_wave,3);
            x1P = sum(abs(x1).^2);

            % Find the index of first receive antenna which has maximum
            % correlation energy
            idx = find(x1P == max(x1P),1);

            % Combine the received waveform which have maximum correlation
            % energy
            x_wave_combined = sum(x1(:,idx(1)),2);

            % Compute FFT of the resultant waveform
            x_fft = fftshift(fft(x_wave_combined));

            % Store the value and location of peak
            [maxValue(shiftIdx),binIndex(shiftIdx)] = max(x_fft);
        end

        % FFT bin values
        fftBinValues = (-rxLen/2:(rxLen/2-1))*(ofdmInfo.SampleRate/rxLen);

        % Find the shift index that corresponds to the maximum of peak
        % value of all the shifted waveforms. Use the FFT bin index
        % corresponding to this maximum shift index. The FFT bin value
        % corresponding to this bin index is the integer frequency offset.
        [~,maxId] = max(maxValue);
        loc = binIndex(maxId);
        out = fftBinValues(loc);
        shiftOut = shiftValues(maxId);
    else
        out = 0;
        shiftOut = 1+sampleOffset;
    end

end

function [out,t] = compensateDopplerShift(inWave,fs,fdSat,flag,t)
% Perform Doppler shift correction

    t1 = (0:size(inWave,1)-1)'/fs;
    if nargin == 5
        t1 = t1 + t; % Add the sample time offset
    end
    if flag
        out = inWave.*exp(1j*2*pi*(-fdSat)*t1);
    else
        out = inWave;
    end
    t = t1(end);

end

Related Topics