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 standards [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

• Doppler pre-compensation at the transmitter, and Doppler compensation at the receiver

• Optional power amplifier modeling with and without memory

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

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 range of transmit power values of the transmit power loop in parallel.

Configure Simulation Length, Transmitter, and Receiver

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 range of transmit power values to simulate. The transmitter power is defined as the power of the time-domain waveform before passing to the power amplifier. The receiver includes its noise figure and the antenna temperature. The noise figure models the receiver internal noise, and the antenna temperature models the input noise. This receiver specifies the noise per antenna element.

simParameters = struct;                   % Clear simParameters variable to contain 
                                          % all key simulation parameters
simParameters.NFrames = 2;                % Number of 10 ms frames
simParameters.TxPower = 60:65;            % Transmit power (dBm)
simParameters.RxNoiseFigure = 6;          % Noise figure (dB)
simParameters.RxAntennaTemperature = 290; % Antenna temperature (K)

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

displaySimulationInformation = true;

Power Amplifier Configuration

The example supports both memory and memoryless power amplifier modeling.

To configure a memory or memoryless power amplifier nonlinearity, use enablePA.

Memoryless Power Amplifier

You can choose one of these memoryless power amplifier models, as defined in Annex A of TR 38.803.

• 2.1 GHz Gallium Arsenide (GaAs)

• 2.1 GHz Gallium Nitride (GaN)

• 28 GHz complementary metal-oxide semiconductor (CMOS)

• 28 GHz GaN

Alternatively, you can set paModel to Custom and use paCharacteristics variable to define the memoryless power amplifier characteristics as a matrix with three columns. The first column defines the input power in dBm. The second column defines the output power in dBm. The third column defines the output phase in degrees. When the paCharacteristics variable is set to empty and the paModel is set to Custom, this example uses a 2.1 GHz laterally-diffused metal-oxide semiconductor (LDMOS) Doherty-based amplifier.

The memoryless nonlinearity applied to the waveform follows this equation for power amplifiers (excluding a custom configuration).

$${\mathit{y}}_{\mathit{P}}\left(\mathit{n}\right)=\sum _{\mathit{k}\in {\mathit{K}}_{\mathit{p}}}{\mathit{a}}_{\mathit{k}}\mathit{x}\left(\mathit{n}\right){\left|\mathit{x}\left(\mathit{n}\right)\right|}^{2\mathit{k}}$$

In this equation,

• ${\mathit{y}}_{\mathit{P}}\left(\mathit{n}\right)$ is the output signal.

• $\mathit{x}\left(\mathit{n}\right)$ is the input signal.

• ${\mathit{K}}_{\mathit{p}}$ is the set of polynomial degree(s).

• ${\mathit{a}}_{\mathit{k}}$ is the polynomial coefficient.

Power Amplifier With Memory

The nonlinearity with memory applied to the waveform follows this equation.

$${\mathit{y}}_{\mathit{P}}\left(\mathit{n}\right)=\sum _{\mathit{m}=0}^{\mathit{M}-1}\sum _{\mathit{k}=0}^{\mathit{K}-1}{\mathit{a}}_{\mathrm{mk}}\mathit{x}\left(\mathit{n}-\mathit{m}\right){\left|\mathit{x}\left(\mathit{n}-\mathit{m}\right)\right|}^{\mathit{k}}$$

In this equation,

• $\mathit{M}$is the memory-polynomial depth.

• $\mathit{K}$is the memory-polynomial degree.

• ${\mathit{a}}_{\mathrm{mk}}$ is the polynomial coefficient.

To model the power amplifier with memory, set the hasMemory variable to true and use the coefficients variable to provide the polynomial coefficients. The coefficients variable is a matrix with number of rows corresponding to memory-polynomial depth and number of columns corresponding to memory-polynomial degree. When coefficients is set to empty, a default value is applied.

By default, the example sets enablePA to false.

enablePA = false;                 % true or false
hasMemory = false;                % true or false
paModel = "2.1GHz GaAs"; % "2.1GHz GaAs", "2.1GHz GaN", "28GHz CMOS", "28GHz GaN", or "Custom"
paCharacteristics = [];         % Lookup table as empty or a matrix with columns: Pin (dBm) | Pout (dBm) | Phase (degrees)
coefficients = [];              % Memory polynomial coefficients

Doppler Compensation Configuration

The example supports two Doppler compensation configurations: one at the transmitter end and the other at the receiver end. For compensation at transmitter end, enable DopplerPreCompensator. Setting the DopplerPreCompensator field to true accounts for Doppler due to satellite movement by applying Doppler pre-compensation to the transmitted waveform. For compensation at receiver end, enable RxDopplerCompensator. Setting the RxDopplerCompensator field to true to estimates and compensates the Doppler shift of 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';

% Include or exclude free space path loss
simParameters.IncludeFreeSpacePathLoss = true;

% 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 transmit power 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 power amplifier nonlinearities — Apply the memory or memoryless power amplifier nonlinearities to the baseband OFDM signal.

5. 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.

6. Model and apply a noisy channel — Pass the generated waveform through an NTN narrowband or NTN TDL fading channel to get the faded waveform. Apply path loss and add thermal noise to the faded waveform.

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

8. 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.

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

10. 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.

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

12. 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.

13. 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.

numTxPowerPoints = length(simParameters.TxPower); % Number of transmit power points
% Array to store the maximum throughput for all transmit power points
maxThroughput = zeros(numTxPowerPoints,1); 
% Array to store the simulation throughput for all transmit power points
simThroughput = zeros(numTxPowerPoints,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;

% Compute the noise amplitude per receive antenna
kBoltz = physconst('Boltzmann');
NF = 10^(simParameters.RxNoiseFigure/10);
T0 = 290;                                               % Noise temperature at the input (K)
Teq = simParameters.RxAntennaTemperature + T0*(NF-1);   % K
N0_ampl = sqrt(kBoltz*waveformInfo.SampleRate*Teq/2.0);

% Compute path loss based on the elevation angle and satellite altitude
re = physconst("earthradius");
c = physconst("lightspeed");
h = simParameters.SatelliteAltitude;
elevAngle = simParameters.ElevationAngle;
d = -re*sind(elevAngle) + sqrt((re*sind(elevAngle)).^2 + h*h + 2*re*h);
lambda = c/simParameters.CarrierFrequency;
pathLoss = fspl(d,lambda)*double(simParameters.IncludeFreeSpacePathLoss); % in dB

% Create a power amplifier System object
hpa = HelperPANonlinearity(HasMemory=hasMemory);
hpaDelay = 0;
if hpa.HasMemory
    if ~isempty(coefficients)
        hpa.Coefficients = coefficients;
    end
    hpaDelay = size(hpa.Coefficients,1)-1;
else
    hpa.Model = paModel;
    if hpa.Model == "Custom" && ~isempty(paCharacteristics)
        hpa.Table = paCharacteristics;
    end
end

% 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 txPowIdx = 1:numTxPowerPoints     % Comment out for parallel computing
% parfor txPowIdx = 1:numTxPowerPoints % 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 transmit power 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 transmit power point
    % experiences the same noise realization
    rng(0,"twister");

    % 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 transmit power point
    pathFilters = [];
    thres = dtxThresold;
    sampleOffset = sampleDelayOffset;
    usePrevShift = usePreviousShift;
    N0 = N0_ampl;
    pl_dB = pathLoss;

    % Transmit power value in dBm
    txPowerdBm = simLocal.TxPower(txPowIdx);

    % 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 transmit power point experiences the
    % same channel realization
    reset(ntnChanLocal.BaseChannel);
    reset(ntnChanLocal.ChannelFilter);

    % Reset the power amplifier
    reset(hpa);

    % 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
        txWaveform0 = nrOFDMModulate(carrier,pdschGrid);

        % Scale the waveform power based on the input transmit power
        wavePower = 10*log10(sum(var(txWaveform0)));
        desiredPower = (txPowerdBm-30)-wavePower;      % In dB
        txWaveform = db2mag(desiredPower)*txWaveform0;

        % Pass the waveform through power amplifier
        if (enablePA == 1)
            txWaveform = hpa(txWaveform);
        end

        % 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);

        % Apply path loss to the signal
        rxWaveform = rxWaveform*db2mag(-pl_dB);

        % Add thermal noise to the received time-domain waveform. Multiply
        % the noise variance with 2 as wgn function performs the scaling
        % within.
        noise = wgn(size(rxWaveform,1),size(rxWaveform,2),2*(N0^2),1,"linear","complex");
        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
        % sum of maximum channel delay and power amplifier delay
        if offset > (maxChDelay+hpaDelay)
            warning(['Estimated timing offset (%d) is greater than the total' ...
                ' delay (%d). This might result in a decoding failure.' ...
                ' The estimated timing offset might be greater than the total' ...
                ' delay because of low transmit power, higher power amplifier distortion,' ...
                ' or not enough DM-RS symbols to synchronize successfully.'], ...
                offset,maxChDelay+hpaDelay);
        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(txPowIdx) = simThroughput(txPowIdx) + sum(~blkerr .* trBlkSizes);
        maxThroughput(txPowIdx) = maxThroughput(txPowIdx) + 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 transmit power %d dBm: %.4f\n',...
            simLocal.NFrames,txPowerdBm,1e-6*simThroughput(txPowIdx)/(simLocal.NFrames*10e-3));
        fprintf('Throughput(%%) for %d frame(s) at transmit power %d dBm: %.4f\n',...
            simLocal.NFrames,txPowerdBm,simThroughput(txPowIdx)*100/maxThroughput(txPowIdx));
    end

end

Throughput(Mbps) for 2 frame(s) at transmit power 60 dBm: 0.0000

Throughput(%) for 2 frame(s) at transmit power 60 dBm: 0.0000

Throughput(Mbps) for 2 frame(s) at transmit power 61 dBm: 0.0000

Throughput(%) for 2 frame(s) at transmit power 61 dBm: 0.0000

Throughput(Mbps) for 2 frame(s) at transmit power 62 dBm: 0.3368

Throughput(%) for 2 frame(s) at transmit power 62 dBm: 5.0000

Throughput(Mbps) for 2 frame(s) at transmit power 63 dBm: 5.7256

Throughput(%) for 2 frame(s) at transmit power 63 dBm: 85.0000

Throughput(Mbps) for 2 frame(s) at transmit power 64 dBm: 6.7360

Throughput(%) for 2 frame(s) at transmit power 64 dBm: 100.0000

Throughput(Mbps) for 2 frame(s) at transmit power 65 dBm: 6.7360

Throughput(%) for 2 frame(s) at transmit power 65 dBm: 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.TxPower,simThroughput*100./maxThroughput,'o-.')
xlabel('Input Transmit Power (dBm)'); 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));

% 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, TxPower = 60:70) 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-compensation is achieved by setting the DopplerPreCompensator field to true and RxDopplerCompensator field to false. The line corresponding to the Rx Doppler Compensation is achieved by setting the DopplerPreCompensator field to false and RxDopplerCompensator field to true.

Further Exploration

You can use this example to further explore these options:

• To analyze the throughput at each transmit power 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:

• HARQEntity — Manage a set of parallel HARQ processes

• HelperGenerateNTNChannel — Generate NTN channel

• HelperPANonlinearity — Model power amplifier nonlinearity

• HelperSetupNTNChannel — Set up NTN channel

• hPRGPrecode — Precoding for PDSCH precoding resource block group (PRG) bundling

• hSkipWeakTimingOffset — Skip timing offset estimates with weak correlation

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

• Model NR NTN Channel

×

Abrir ejemplo

Tiene una versión modificada de este ejemplo. ¿Desea abrir este ejemplo con sus modificaciones?

No, sobrescribir la versión modificada Sí