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 hybrid automatic repeat request (HARQ) support up to 32 processes.

• Optional power amplifier modeling with and without memory. The power amplifier modeling with memory requires RF Toolbox™.

• Optional static and time-varying propagation delay modeling.

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 performing Doppler pre-compensation and includes the gain of 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;                   % Create simParameters structure 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 model the power amplifier with memory, this example requires RF Toolbox™.

To configure a memory or memoryless power amplifier nonlinearity, use the enablePA variable. The input signal that is passed to the power amplifier is a normalized signal with the maximum signal amplitude.

Memoryless Power Amplifier

You can select one of these memoryless power amplifier models (paModel), 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 you set the paCharacteristics variable to empty and the paModel to Custom, this example uses a 2.1 GHz laterally-diffused metal-oxide semiconductor (LDMOS) Doherty-based amplifier.

When you set paModel to a value other than Custom, the memoryless nonlinearity applied to the waveform follows this equation for power amplifiers.

$${\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 memory polynomial 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 you set coefficients 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

When you set enablePA is true, use the scaleFactor variable to modify the maximum input signal amplitude to excite the power amplifier nonlinearity. scaleFactor controls the operating region of the power amplifier and is applied for each transmit antenna. You can also use the scaleFactor variable to set power backoff. For example, to provide a power backoff of 3 dB to a signal passed through the power amplifier, set scaleFactor to -3. Ensure the input signal is within the characterization range of the power amplifier model.

When scaleFactor is empty, the example uses a default value of -35 dB in these cases.

• hasMemory is false, paModel is Custom, and paCharacteristics is empty.

• hasMemory is true and coefficients is empty.

In all other cases, when you set scaleFactor to empty, the example uses a default value of 0 dB.

scaleFactor = []; % Amplitude scaling, in dB

Doppler Compensation Configuration

The example supports two Doppler compensation configurations: one at the transmitter and the other at the receiver. For compensation at the transmitter, 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 the receiver, enable the RxDopplerCompensator field. Setting the RxDopplerCompensator field to true estimates and compensates the Doppler shift of the received waveform, using cyclic prefix and reference signals. When you set the RxDopplerCompensator field to true, you can select the technique for estimating and compensating Doppler at the receiver using the RxDopplerCompensationMethod field. The RxDopplerCompensationMethod field supports:

• Independent time-frequency synchronization (independent time-freq), where the receiver compensates for frequency or Doppler shift first and then compensates for timing offset.

• Joint time-frequency synchronization (joint time-freq), where the receiver compensates for both frequency and time at once.

By default, this example assumes the user equipment (UE) is at the satellite beam center. To model a UE in the satellite beam other than the beam center (as shown in the next image) and to apply a Doppler shift common to all the UEs in a beam (${\mathit{f}}_{\mathit{d},\mathrm{common}}$), you can use the PreCompensationDopplerShift field. When the PreCompensationDopplerShift field is empty, the example uses the Doppler shift due to satellite at the UE (${\mathit{f}}_{\mathit{d},\mathrm{sat}}$) as ${\mathit{f}}_{\mathit{d},\mathrm{common}}$. This example assumes ${\mathit{f}}_{\mathit{d},\mathrm{common}}$ is known. To observe the link performance when PreCompensationDopplerShift field is nonempty, both DopplerPreCompensator field and RxDopplerCompensator field must be set to true. Enabling DopplerPreCompensator compensates ${\mathit{f}}_{\mathit{d},\mathrm{common}}$, while enabling RxDopplerCompensator compensates the residual satellite Doppler shift (${\mathit{f}}_{\mathit{d},\mathrm{sat}}$-${\mathit{f}}_{\mathit{d},\mathrm{common}}$) along with Doppler shift due to UE movement.

simParameters.DopplerPreCompensator = true;
simParameters.PreCompensationDopplerShift = [];                   % In Hz
simParameters.RxDopplerCompensator = false;
simParameters.RxDopplerCompensationMethod = "independent time-freq";
% The example uses below fields to estimate Doppler shift, when
% RxDopplerCompensator is set to true and RxDopplerCompensationMethod is
% set to joint time-freq.
% Set the search range of Doppler shift in Hz [MIN,MAX]
simParameters.FrequencyRange = [-50e3 50e3];
% Set the search range resolution of Doppler shift in Hz
simParameters.FrequencyResolution = 1e3;

Initial Timing Synchronization Algorithm Selection

Select an algorithm for initial timing synchronization.

• Auto correlation (auto corr): The receiver performs timing synchronization using auto correlation with PDSCH DM-RS.

• Differential correlation (diff corr): The receiver performs timing synchronization using differential correlation with PDSCH DM-RS.

• Joint time-frequency technique (joint time-freq): The receiver performs timing synchronization by compensating for initial frequency and initial timing at once.

simParameters.InitialTimingSynchronization = "joint time-freq";
% The example uses below fields to perform initial synchronization, when
% InitialTimingSynchronization is set to joint time-freq.
% Set the initial search range of Doppler shift in Hz [MIN,MAX]
simParameters.InitialFrequencyRange = [-50e3 50e3];
% Set the initial search range resolution of Doppler shift in Hz
simParameters.InitialFrequencyResolution = 1e3;

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, 120, 240, 480, or 960

• Cyclic prefix length (CP): normal or extended

• Cell identity

• Number of transmit and receive antennas

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 = []; % Use empty to auto-configure the DM-RS ports
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;
% 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
% Define data type for resource grids and waveforms
simParameters.DataType = "double";

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 an NTN TDL channel.
simParameters.NTNChannelType = "Narrowband";

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

% Delay model configuration
% This example models only one-way propagation delay and provides immediate
% feedback without any delay
simParameters.DelayModel = "None";    % "None", "Static", or "Time-varying"

% 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.MobileAltitude = 0;                      % Mobile altitude (in m)
simParameters.SatelliteAltitude = 600000;              % Satellite altitude (in m)
simParameters.SampleRate = waveformInfo.SampleRate;
simParameters.RandomStream = "mt19937ar with seed";
simParameters.Seed = 73;
simParameters.OutputDataType = simParameters.DataType;

% Set the following fields for NTN narrowband channel
if simParameters.NTNChannelType == "Narrowband"
    simParameters.Environment = "Urban";
    simParameters.AzimuthOrientation = 0;
end

% Set the following fields for NTN TDL channel
if simParameters.NTNChannelType == "TDL"
    simParameters.DelayProfile = "NTN-TDL-A";
    simParameters.DelaySpread = 30e-9;
end

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

% Calculate the Doppler shift due to satellite movement
c = physconst("lightspeed");
satelliteDopplerShift = dopplerShiftCircularOrbit( ...
    simParameters.ElevationAngle,simParameters.SatelliteAltitude, ...
    simParameters.MobileAltitude,simParameters.CarrierFrequency);

% Define NTN narrowband channel based on the specified fields in
% simParameters structure
if simParameters.NTNChannelType == "Narrowband"
    channel = p681LMSChannel;
    channel.Environment = simParameters.Environment;
    channel.AzimuthOrientation = simParameters.AzimuthOrientation;
    channel.CarrierFrequency = simParameters.CarrierFrequency;
    channel.ElevationAngle = simParameters.ElevationAngle;
    channel.MobileSpeed = simParameters.MobileSpeed;
    channel.SatelliteDopplerShift = satelliteDopplerShift;
end

% Define NTN TDL channel based on specified fields in simParameters
% structure
if simParameters.NTNChannelType == "TDL"
    channel = nrTDLChannel;
    channel.DelayProfile = simParameters.DelayProfile;
    channel.DelaySpread = simParameters.DelaySpread;
    channel.SatelliteDopplerShift = satelliteDopplerShift;
    channel.MaximumDopplerShift = ...
        simParameters.MobileSpeed*simParameters.CarrierFrequency/c;
    channel.NumTransmitAntennas = simParameters.NumTransmitAntennas;
    channel.NumReceiveAntennas = simParameters.NumReceiveAntennas;
end

% Assign the parameters common to both TDL and narrowband channels
channel.SampleRate = simParameters.SampleRate;
channel.RandomStream = simParameters.RandomStream;
channel.Seed = simParameters.Seed;

% 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 buffer
% and process the received signal with the expected length.
chInfo = info(channel);
maxChDelay = ceil(max(chInfo.PathDelays*channel.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 provides the time-domain waveform by performing OFDM modulation of the generated resource grid. Normalize the waveform with the maximum waveform amplitude for each antenna.

4. Apply power amplifier nonlinearities — Scale the amplitude of normalized waveform depending on the input scaling factor. Apply the memory or memoryless power amplifier nonlinearities to the baseband OFDM waveform. Scale the waveform power to the desired transmit power.

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 — Delay the generated waveform depending on the propagation latency. Pass the delayed 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. Perform initial synchronization — Check the presence of signal using energy detection. After energy detection, the received waveform uses the PDSCH DM-RS of initial slot to get the initial timing offset. Perform this step till the initial synchronization is achieved.

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

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

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

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

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

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

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

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

% Number of transmit power points
numTxPowerPoints = length(simParameters.TxPower);
% 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);
% Array to store the signal-to-noise ratio (SNR) for all transmit power points
snrVec = zeros(numTxPowerPoints,1);

% Common Doppler shift for use in the simulations
if isempty(simParameters.PreCompensationDopplerShift)
    commonDopplerShift = satelliteDopplerShift;
else
    commonDopplerShift = simParameters.PreCompensationDopplerShift;
end

% 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 path loss based on the elevation angle and satellite altitude
slantDist = slantRangeCircularOrbit(simParameters.ElevationAngle, ...
    simParameters.SatelliteAltitude,simParameters.MobileAltitude);
lambda = c/simParameters.CarrierFrequency;
pathLoss = fspl(slantDist,lambda)*double(simParameters.IncludeFreeSpacePathLoss); % in dB

% Get the slot time based on subcarrier spacing
symLen = cumsum(waveformInfo.SymbolLengths);
nSymbSlot = simParameters.Carrier.SymbolsPerSlot;
nSubFrames = 10;
samples = [0 symLen(nSymbSlot:nSymbSlot:end-1)]'./waveformInfo.SampleRate;
subframeTimes = (0:(simParameters.NFrames*nSubFrames)-1)*1e-3;
slotTimes = cast(reshape(samples+subframeTimes,[],1),simParameters.DataType);

% Maximum variable propagation delay
maxVarPropDelay = 0;

% Compute static delay in seconds and samples
delayInSeconds = cast(slantDist./c,simParameters.DataType);
delayInSamples = delayInSeconds.*waveformInfo.SampleRate;
integDelaySamples = floor(delayInSamples);
fracDelaySamples = (delayInSamples - integDelaySamples);
numVariableFracDelaySamples = repmat(fracDelaySamples,1, ...
    (simParameters.Carrier.SlotsPerFrame*simParameters.NFrames) + 1);
numVariableIntegSamples = 0;
pathLoss = repmat(pathLoss,1,numel(slotTimes));
% Initialize configuration objects for delay modeling
if simParameters.DelayModel == "None"
    staticDelay = dsp.Delay(Length=0);
    variableIntegerDelay = 0;
    variableFractionalDelay = 0;
    delayInSeconds = 0;
elseif simParameters.DelayModel == "Static"
    staticDelay = dsp.Delay(Length=integDelaySamples);
    variableIntegerDelay = 0;
    variableFractionalDelay = dsp.VariableFractionalDelay(...
        InterpolationMethod="Farrow", ...
        FarrowSmallDelayAction="Use off-centered kernel", ...
        MaximumDelay=1);
else
    % Model time-varying delay assuming satellite moves in a circular orbit.
    % The example applies delay for each slot and assumes there is no
    % significant change in delay for each sample.

    % Calculate delay across the simulation time in steps of slot time
    SU = slantRangeCircularOrbit(simParameters.ElevationAngle, ...
        simParameters.SatelliteAltitude,simParameters.MobileAltitude,slotTimes);
    pathLoss = fspl(SU,lambda).*double(simParameters.IncludeFreeSpacePathLoss);
    delayInSeconds = SU./c;
    delayInSamples = delayInSeconds.*waveformInfo.SampleRate;

    % Compute dynamic range of delay and configure the delay objects
    % accordingly
    integDelaySamples = floor(delayInSamples);
    numStaticDelaySamples = min(integDelaySamples);
    remVariableDelaySamples = delayInSamples - numStaticDelaySamples;
    staticDelay = dsp.Delay(Length=numStaticDelaySamples);
    numVariableIntegSamples = floor(remVariableDelaySamples);
    numVariableIntegSamples(numVariableIntegSamples < 0) = 0;
    maxVarPropDelay = max(numVariableIntegSamples)+2;
    variableIntegerDelay = dsp.VariableIntegerDelay(...
        MaximumDelay=maxVarPropDelay);
    numVariableFracDelaySamples = remVariableDelaySamples - numVariableIntegSamples;
    variableFractionalDelay = dsp.VariableFractionalDelay( ...
        InterpolationMethod="Farrow", ...
        FarrowSmallDelayAction="Use off-centered kernel", ...
        MaximumDelay=1);
end

% Check the number of HARQ processes and initial propagation delay
initialSlotDelay = find(slotTimes>=delayInSeconds(1),1)-1;
if simParameters.PDSCHExtension.EnableHARQ
    if simParameters.PDSCHExtension.NHARQProcesses < initialSlotDelay
        error("In case of HARQ, this example supports transmission of continuous data only. " + ... 
            "Set the number of HARQ processes (" + (simParameters.PDSCHExtension.NHARQProcesses) +...
            ") to a value greater than or equal to the maximum propagation delay in slots (" + ...
            initialSlotDelay +").")
    end
end

% Initial frequency shift search space
if simParameters.InitialTimingSynchronization == "joint time-freq"
    inifVals = simParameters.InitialFrequencyRange(1):simParameters.InitialFrequencyResolution:simParameters.InitialFrequencyRange(2);
else
    inifVals = 0;
end

% Frequency shift search space
if simParameters.RxDopplerCompensator == 1 ...
        && simParameters.RxDopplerCompensationMethod == "joint time-freq"
    fVals = simParameters.FrequencyRange(1):simParameters.FrequencyResolution:simParameters.FrequencyRange(2);
else
    % In case of no receiver Doppler compensation, treat the frequency
    % value is 0 Hz.
    fVals = 0;
end

% Initialize the power amplifier function handle or System object depending
% on the input configuration
paInputScaleFactor = 0;                                 % in dB
hpaDelay = 0;
if hasMemory == 1
    hpa = rf.PAmemory;                                  % Requires RF Toolbox
    if isempty(coefficients)
        hpa.CoefficientMatrix = getDefaultCoefficients;
        paInputScaleFactor = -35;
    else
        hpa.CoefficientMatrix = coefficients;
    end
    hpa.UnitDelay = 1;
    hpaDelay = size(hpa.CoefficientMatrix,1)-1;
else
    if paModel == "Custom"
        if isempty(paCharacteristics)
            hpa = comm.MemorylessNonlinearity(Method="Lookup table", ...
                Table=getDefaultLookup);
            paInputScaleFactor = -35;
        else
            hpa = comm.MemorylessNonlinearity(Method="Lookup table", ...
                Table=paCharacteristics);
        end
    elseif paModel == "2.1GHz GaAs"
        hpa = @(in) paMemorylessGaAs2Dot1GHz(in);
    elseif paModel == "2.1GHz GaN"
        hpa = @(in) paMemorylessGaN2Dot1GHz(in);
    elseif paModel == "28GHz CMOS"
        hpa = @(in) paMemorylessCMOS28GHz(in);
    else % "28GHz GaN"
        hpa = @(in) paMemorylessGaN28GHz(in);
    end
end
% Repeat hpa to have independent processing for each antenna
hpa = repmat({hpa},1,simParameters.NumTransmitAntennas);

% Update the power amplifier input scaling factor, based on scaleFactor
if ~isempty(scaleFactor)
    paInputScaleFactor = scaleFactor;
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 usePreviousShift variable 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;

% Set useDiffCorr variable to true, to use the shift estimated from
% differential correlation directly in the integer Doppler shift
% estimation. When set to false, the range of shift values also include the
% shift estimated from differential correlation.
useDiffCorr = true;

% Set the amplitude scaling factor to use in energy detection. For the
% default case, a factor of 1.03 is used to avoid missed detections at 60
% dBm transmit power. A value of 0 assumes each sample is an actual signal.
amplThreshold = 1.03;

% Use the minimum number of samples for a slot in the whole frame as
% window length
mrms = dsp.MovingRMS;
slotsPerSubFrameFlag = simParameters.Carrier.SlotsPerSubframe > 1;
mrms.WindowLength = symLen((1+(slotsPerSubFrameFlag))*simParameters.Carrier.SymbolsPerSlot) ...
    -slotsPerSubFrameFlag*symLen(simParameters.Carrier.SymbolsPerSlot);

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

    % Make copies of channel-level parameters to simplify subsequent
    % parameter referencing
    carrier = simLocal.Carrier;
    rxCarrier = 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

    % Make copies of intermediate variables to have warning-free execution
    % with Parallel Computing Toolbox
    thres = dtxThresold;
    sampleOffset = sampleDelayOffset;
    usePrevShift = usePreviousShift;
    useDiffCorrFlag = useDiffCorr;
    N0 = N0_ampl;
    pl_dB = pathLoss;
    varIntegSamples = numVariableIntegSamples;
    varFracSamples = numVariableFracDelaySamples;
    fValsVec = fVals;
    inifValsVec = inifVals;
    threshFactor = amplThreshold;
    initialDelay = initialSlotDelay;
    preDopplerShift = commonDopplerShift;

    % Initialize temporary variables
    offset = 0;
    shiftOut = 0;
    txHarqProc = 0;
    rxHarqProc = 0;
    prevWave = [];
    pathFilters = [];
    rxBuff = [];
    syncCheck = true;

    % Reset the channel so that each transmit power point experiences the
    % same channel realization
    reset(channel);

    % Reset the power amplifier
    for numHPA = 1:numel(hpa)
        if ~isa(hpa{numHPA},"function_handle")
            reset(hpa{numHPA})
        end
    end

    % Reset the delay objects
    reset(staticDelay)
    if isa(variableIntegerDelay,"dsp.VariableIntegerDelay")
        reset(variableIntegerDelay)
    end
    if isa(variableFractionalDelay,"dsp.VariableFractionalDelay")
        reset(variableFractionalDelay)
    end

    % Reset the moving RMS object
    reset(mrms)

    % 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
    % Create a parallel array of all HARQ processes
    harqEntity = cell(pdschextra.NHARQProcesses,1);
    for harqId = 1:pdschextra.NHARQProcesses
        harqEntity{harqId} = HARQEntity(harqSequence(harqId),rvSeq,pdsch.NumCodewords);
    end

    % 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,channel,simLocal.DataType);
    newWtx = getPrecodingMatrix(carrier,pdsch,estChannelGrid,pdschextra.PRGBundleSize);

    % 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{txHarqProc+1}.NewData(cwIdx)
                trBlk = randi([0 1],trBlkSizes(cwIdx),1,'int8');
                setTransportBlock(encodeDLSCH,trBlk,cwIdx-1,harqEntity{txHarqProc+1}.HARQProcessID);
                % Flush decoder soft buffer explicitly for any new data
                % because of previous RV sequence time out
                if harqEntity{txHarqProc+1}.SequenceTimeout(cwIdx)
                    resetSoftBuffer(decodeDLSCHLocal,cwIdx-1,harqEntity{txHarqProc+1}.HARQProcessID);
                end
            end
        end

        % Encode the DL-SCH transport blocks
        codedTrBlocks = encodeDLSCH(pdsch.Modulation,pdsch.NumLayers, ...
            pdschIndicesInfo.G,harqEntity{txHarqProc+1}.RedundancyVersion, ...
            harqEntity{txHarqProc+1}.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, ...
            OutputDataType=simLocal.DataType);

        % Perform implementation-specific PDSCH MIMO precoding and mapping
        [pdschAntSymbols,pdschAntIndices] = nrPDSCHPrecode( ...
            carrier,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] = nrPDSCHPrecode( ...
            carrier,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] = nrPDSCHPrecode( ...
            carrier,ptrsSymbols,ptrsIndices,wtx);
        pdschGrid(ptrsAntIndices) = ptrsAntSymbols;

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

        % Normalize the waveform with maximum waveform amplitude
        txWaveform = txWaveform0./max(abs(txWaveform0));

        % Adjust the waveform amplitude and pass the waveform through power
        % amplifier
        if (enablePA == 1)
            % Scale the amplitude of the waveform, as applicable
            txWaveform = txWaveform.*db2mag(paInputScaleFactor);

            % Pass the adjusted waveform through the power amplifier
            for colIdx = 1:size(txWaveform,2)
                hpaTemp = hpa{colIdx};
                txWaveform(:,colIdx) = hpaTemp(txWaveform(:,colIdx));
            end
        end

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

        % Apply Doppler pre-compensation depending on DopplerPreCompensator
        % field
        txWaveform = compensateDopplerShift(...
            txWaveform,channel.SampleRate, ...
            preDopplerShift,simLocal.DopplerPreCompensator);

        % Apply path loss to the signal
        txWaveform = txWaveform*db2mag(-pl_dB(carrier.NSlot+1));

        % Apply fixed or static delay
        delayedTx = staticDelay(txWaveform);
        % Apply variable integer delay
        if isa(variableIntegerDelay,"dsp.VariableIntegerDelay")
            delayedTx = variableIntegerDelay(delayedTx,varIntegSamples(carrier.NSlot+1));
        end
        % Apply variable fractional delay
        if isa(variableFractionalDelay,"dsp.VariableFractionalDelay")
            delayedTx = variableFractionalDelay(delayedTx,varFracSamples(carrier.NSlot+1));
        end

        % Pass the waveform through the channel
        txWaveform = delayedTx;
        [rxWaveform,pathGains] = channel(txWaveform);

        % 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");
        sigPowerRE = sum(var(rxWaveform))*((waveinfoLocal.Nfft)^2/(carrier.NSizeGrid*12));
        noisePowerRE = sum(var(noise))*(waveinfoLocal.Nfft);
        snrVec(txPowIdx) = snrVec(txPowIdx) + sigPowerRE./noisePowerRE;
        rxWaveform = rxWaveform + cast(noise,simLocal.DataType);

        % Update the transmit HARQ process number
        if pdschextra.EnableHARQ
            txHarqProc = mod(txHarqProc+1,pdschextra.NHARQProcesses);
        end

        % Compute the moving RMS of the signal and perform energy detection
        % for initial synchronization
        metric = mrms(complex(rxWaveform));
        idx = metric > (sqrt(2)*N0*threshFactor);
        if ~any(idx(:)) && (rxCarrier.NSlot == 0)
            % Store the waveform that didn't pass the metric to use for
            % initial synchronization
            prevWave = rxWaveform;
            continue;
        end
        % Provide a warning when initial synchronization is missed
        if (rxCarrier.NSlot == 0) && syncCheck
            syncCheck = false;
            if (carrier.NSlot > initialDelay)
                warning("Initial slot synchronization is missed for transmit power of %d dBm. " + ...
                    "This can cause failure of all the slots. " + ...
                    "For proper synchronization, increase the transmit power.",txPowerdBm)
            end
        end

        % Buffer all the valid signal such that the length of 3 slots is
        % used for initial synchronization, and some portion of previous
        % slot is used for next slot.
        rxBuff = [rxBuff;prevWave;rxWaveform]; %#ok<AGROW>
        prevWave = [];
        if (size(rxBuff,1) < (3*(mrms.WindowLength))) && (rxCarrier.NSlot == 0)
            continue
        else
            % Here onwards reception happens continuously

            % Use the whole buffered waveform for receiver and generate the
            % reference signals for this particular slot to use for
            % waveform processing.
            rxData = rxBuff;
            refDMRSSymbols = nrPDSCHDMRS(rxCarrier,pdsch);
            refDMRSIndices = nrPDSCHDMRSIndices(rxCarrier,pdsch);
            refPTRSSymbols = nrPDSCHPTRS(rxCarrier,pdsch);
            refPTRSIndices = nrPDSCHPTRSIndices(rxCarrier,pdsch);
        end

        % Gather the number of samples to be processed in current
        % receiver slot.
        % 1. Find the number of cyclic prefix samples used in the current
        % receiver slot and sum of all samples
        % 2. Add the FFT size corresponding to number of OFDM symbols
        % in current slot with the resultant value in step 1
        cpl = circshift(waveinfoLocal.CyclicPrefixLengths,-rxCarrier.NSlot*rxCarrier.SymbolsPerSlot);
        numSamplesInRxSlot = waveinfoLocal.Nfft*rxCarrier.SymbolsPerSlot + sum(cpl(1:rxCarrier.SymbolsPerSlot));

        % Due to large Doppler shift, the estimate using the DM-RS
        % correlation gives an inaccurate estimate. Thus, perform joint
        % time and Doppler shift estimation for the received signal to get
        % the initial timing offset.
        if rxCarrier.NSlot == 0
            if simLocal.InitialTimingSynchronization == "auto corr"
                initialOffset = nrTimingEstimate(rxCarrier,rxData,refDMRSIndices,refDMRSSymbols);
            elseif simLocal.InitialTimingSynchronization == "diff corr"
                initialOffset = diffcorr(rxCarrier,rxData,refDMRSIndices,refDMRSSymbols);
            else
                initialOffset = jointTimeFreq(rxCarrier,rxData,refDMRSIndices,refDMRSSymbols,...
                    inifValsVec);
            end
            d = maxVarPropDelay;
            if d > initialOffset
                d = initialOffset;
            end
        else
            % Use the timing estimate for the required number of samples
            initialOffset = 0;
            d = 0;
        end

        % From the starting position provided by initial offset, consider
        % the length of received waveform such that all the delays due to
        % channel, power amplifier, and propagation distance are covered.
        totalDelay = maxChDelay+maxVarPropDelay+hpaDelay;
        endIdx = initialOffset+numSamplesInRxSlot+totalDelay;
        if endIdx > size(rxData,1)
            endIdx = size(rxData,1);
        end
        rxWaveform = rxData(initialOffset+1:endIdx,:);
        % Update the buffer with portion of present slot data to process
        % the next slot
        rxBuff = rxData(initialOffset-d+(numSamplesInRxSlot+1):end,:);

        if simLocal.RxDopplerCompensator && ...
                simLocal.RxDopplerCompensationMethod == "joint time-freq"
            % Perform joint time-frequency synchronization
            [offset,fOEst] = jointTimeFreq(rxCarrier,rxWaveform,refDMRSIndices,refDMRSSymbols,...
                    fValsVec);
            % Compensate Doppler shift
            rxWaveform = compensateDopplerShift(rxWaveform,waveinfoLocal.SampleRate, ...
                fOEst,true);
            % Estimate and compensate the residual Doppler shift
            [fractionalDopplerShift,detFlag] = estimateFractionalDopplerShift( ...
                rxWaveform,rxCarrier.SubcarrierSpacing,waveinfoLocal.Nfft, ...
                waveinfoLocal.CyclicPrefixLengths(2),0,true);
            rxWaveform = compensateDopplerShift(rxWaveform,waveinfoLocal.SampleRate, ...
                fractionalDopplerShift,true);
            % Get the estimated Doppler shift value
            estimatedDS = fractionalDopplerShift + fOEst;
        else
            % 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,rxCarrier.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( ...
                rxCarrier,rxWaveform,refDMRSIndices,refDMRSSymbols,sampleOffset, ...
                usePrevShift,useDiffCorrFlag,shiftOut-sampleOffset,totalDelay, ...
                (simLocal.RxDopplerCompensator && detFlag));
            rxWaveform = compensateDopplerShift(rxWaveform,waveinfoLocal.SampleRate, ...
                integerDopplerShift,simLocal.RxDopplerCompensator);

            % Get the estimated Doppler shift value
            estimatedDS = fractionalDopplerShift + integerDopplerShift;

            % 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(rxCarrier,rxWaveform, ...
                refDMRSIndices,refDMRSSymbols);
            offset = hSkipWeakTimingOffset(offset,t,mag);
        end
        rxWaveform = rxWaveform(1+offset:end,:);
        if size(rxWaveform,1) > numSamplesInRxSlot
            rxWaveform = rxWaveform(1:numSamplesInRxSlot,:);
        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(rxCarrier,rxWaveform);
        [K,L,R] = size(rxGrid);
        if (L < rxCarrier.SymbolsPerSlot)
            rxGrid = cat(2,rxGrid,zeros(K,rxCarrier.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(rxCarrier,rxGrid,...
            refDMRSIndices,refDMRSSymbols,'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(...
            rxCarrier,estChannelGrid,conj(wtx));

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

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

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

            % Extract PT-RS symbols from received grid and estimated
            % channel grid
            [ptrsRx,ptrsHest,~,~,~,ptrsLayerIndices] = ...
                nrExtractResources(refPTRSIndices,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,refPTRSIndices,refPTRSSymbols,...
                CyclicPrefix=rxCarrier.CyclicPrefix);

            % 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

        if ~any(isnan(pdschEq))
            % Decode PDSCH symbols
            [dlschLLRs,rxSymbols] = nrPDSCHDecode(rxCarrier,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{rxHarqProc+1}.RedundancyVersion,...
                harqEntity{rxHarqProc+1}.HARQProcessID);
        else
            blkerr = true;
        end

        % 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 increment slot number
        updateProcess(harqEntity{rxHarqProc+1},blkerr,trBlkSizes,pdschIndicesInfo.G);
        rxCarrier.NSlot = rxCarrier.NSlot + 1;

        % Increment the receiver HARQ process number
        if pdschextra.EnableHARQ
            rxHarqProc = mod(rxHarqProc+1,pdschextra.NHARQProcesses);
        end

    end

    % Display the results
    if displaySimulationInformation == 1
        snrdb = pow2db(snrVec(txPowIdx)./NSlots);
        if rxCarrier.NSlot == 0
            fprintf("\nNo slot is processed in the receiver at transmit power %d dBm (modeled SNR per RE %.4f dB)", ...
                txPowerdBm,snrdb);
        else
            fprintf("\nThroughput(Mbps) for %d frame(s) at transmit power %d dBm (modeled SNR per RE %.4f dB): %.4f\n",...
                simLocal.NFrames,txPowerdBm,snrdb,1e-6*simThroughput(txPowIdx)/(simLocal.NFrames*10e-3));
            fprintf("Throughput(%%) for %d frame(s) at transmit power %d dBm (modeled SNR per RE %.4f dB): %.4f\n",...
                simLocal.NFrames,txPowerdBm,snrdb,simThroughput(txPowIdx)*100/maxThroughput(txPowIdx));
        end
    end

end

Throughput(Mbps) for 2 frame(s) at transmit power 60 dBm (modeled SNR per RE 5.3762 dB): 0.0000

Throughput(%) for 2 frame(s) at transmit power 60 dBm (modeled SNR per RE 5.3762 dB): 0.0000

Throughput(Mbps) for 2 frame(s) at transmit power 61 dBm (modeled SNR per RE 6.3762 dB): 0.0000

Throughput(%) for 2 frame(s) at transmit power 61 dBm (modeled SNR per RE 6.3762 dB): 0.0000

Throughput(Mbps) for 2 frame(s) at transmit power 62 dBm (modeled SNR per RE 7.3762 dB): 0.3368

Throughput(%) for 2 frame(s) at transmit power 62 dBm (modeled SNR per RE 7.3762 dB): 5.2632

Throughput(Mbps) for 2 frame(s) at transmit power 63 dBm (modeled SNR per RE 8.3762 dB): 5.8940

Throughput(%) for 2 frame(s) at transmit power 63 dBm (modeled SNR per RE 8.3762 dB): 92.1053

Throughput(Mbps) for 2 frame(s) at transmit power 64 dBm (modeled SNR per RE 9.3762 dB): 6.3992

Throughput(%) for 2 frame(s) at transmit power 64 dBm (modeled SNR per RE 9.3762 dB): 100.0000

Throughput(Mbps) for 2 frame(s) at transmit power 65 dBm (modeled SNR per RE 10.3762 dB): 6.3992

Throughput(%) for 2 frame(s) at transmit power 65 dBm (modeled SNR per RE 10.3762 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.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, PreCompensationDopplerShift to [], and RxDopplerCompensator field to false. The line corresponding to the Rx Doppler Compensation is achieved by setting the DopplerPreCompensator field to false, RxDopplerCompensator field to true, and RxDopplerCompensationMethod field to "independent time-freq".

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. Select either "joint time-freq" or "independent time-freq" method for receiver Doppler compensation.

• 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 check the throughput performance in the presence of propagation delay, set the DelayModel to "Static" or "Time-varying".

• 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).

Supporting Files

The example uses these helper functions:

• HARQEntity — Manage a set of parallel HARQ processes

• 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,channel,dataType)
% Obtain channel estimate before first transmission. Use this function to
% obtain a precoding matrix for the first slot.

    ofdmInfo = nrOFDMInfo(carrier);

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

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

    % Filter through channel and get the path gains and path filters
    [~,pathGains,sampleTimes] = channel(tmpWaveform);
    if isa(channel,'nrTDLChannel')
        pathFilters = getPathFilters(channel);
    else
        pathFilters = chInfo.ChannelFilterCoefficients.';
    end

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

    % Perfect channel estimate
    estChannelGrid = nrPerfectChannelEstimate(...
        carrier,pathGains,pathFilters,offset,double(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],'like',hestGrid);
    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 = nrPDSCHPrecode( ...
        carrier,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,'like',rxWave)];

    % 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,useDiffCorr,shiftIn,maxOffset,flag)
% Estimate the integer Doppler shift using demodulation reference signal

    if flag

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

        % Find the timing offset using differential correlation
        offset = diffcorr(carrier,rx,refInd,refSym);
        if offset > maxOffset
            offset = 0;
        end

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

        % Initialize temporary variables
        shiftLen = length(shiftValues);
        maxValue = complex(zeros(shiftLen,1));
        binIndex = zeros(shiftLen,1);
        [rxLen,R] = size(rx);
        xWave = zeros([rxLen P],'like',rx);

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

        % 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 = rx(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
                    xWave(:,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(xWave,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
            xWaveCombined = sum(x1(:,idx(1)),2);

            % Compute FFT of the resultant waveform
            xWaveCombinedTemp = buffer(xWaveCombined,ofdmInfo.Nfft);
            xFFT = sum(fftshift(fft(xWaveCombinedTemp)),2);

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

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

        % 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 = compensateDopplerShift(inWave,fs,fdSat,flag)
% Perform Doppler shift correction

    t = (0:size(inWave,1)-1)'/fs;
    if flag
        out = inWave.*exp(1j*2*pi*(-fdSat)*t);
    else
        out = inWave;
    end

end

function [offset,mag] = diffcorr(carrier,rx,refInd,refSym)
% Perform differential correlation for the received signal

    % Get the number of subcarriers, OFDM symbols, and layers 
    K = carrier.NSizeGrid*12;                % Number of subcarriers
    L = carrier.SymbolsPerSlot;              % Number of OFDM symbols in slot
    P = ceil(max(double(refInd(:))/(K*L)));  % Number of layers

    % Generate the reference signal
    refGrid = nrResourceGrid(carrier,P);
    refGrid(refInd) = refSym;
    refWave = nrOFDMModulate(carrier,refGrid,'Windowing',0);

    % Get the differential of the received signal and reference signal
    waveform = conj(rx(1:end-1,:)).*rx(2:end,:);
    ref = conj(refWave(1:end-1,:)).*refWave(2:end,:);
    [T,R] = size(waveform);

    % To normalize the xcorr behavior, pad the input waveform to make it
    % longer than the reference signal
    refLen = size(ref,1);
    waveformPad = [waveform; zeros([refLen-T R],'like',waveform)];

    % Store correlation magnitude for each time sample, receive antenna and
    % port
    mag = zeros([max(T,refLen) R P],'like',waveformPad);
    for r = 1:R
        for p = 1:P
            % Correlate the given antenna of the received signal with the
            % given port of the reference signal
            refcorr = xcorr(waveformPad(:,r),ref(:,p));
            mag(:,r,p) = abs(refcorr(T:end));
        end
    end

    % Sum the magnitudes of the ports
    mag = sum(mag,3);

    % Find timing peak in the sum of the magnitudes of the receive antennas
    [~,peakindex] = max(sum(mag,2));
    offset = peakindex - 1;

end

function [tO,fO] = jointTimeFreq(carrier,rx,refInd,refSym,fSearchSpace)
% Perform joint time-frequency synchronization
    numFreqVals = length(fSearchSpace);
    peakVal = zeros(numFreqVals,1);
    peakIdx = peakVal;
    ofdmInfo = nrOFDMInfo(carrier);
    fs = ofdmInfo.SampleRate;
    for fIdx = 1:numFreqVals
        rxCorrected = compensateDopplerShift(rx,fs,fSearchSpace(fIdx),true);
        [~,corr] = nrTimingEstimate(carrier,rxCorrected,refInd,refSym);
        corr = sum(abs(corr),2);
        [peakVal(fIdx),peakIdx(fIdx)] = max(corr);
    end
    [~,id] = max(peakVal);
    % Estimate frequency shift and timing offset
    fO = fSearchSpace(id);
    tO = peakIdx(id)-1;
end

% Functions to model power amplifier nonlinearity
function out = paMemorylessGaAs2Dot1GHz(in)
    % 2.1 GHz GaAs
    absIn = abs(in).^(2*(1:7));
    out = (-0.618347-0.785905i) * in + (2.0831-1.69506i) * in .* absIn(:,1) + ...
        (-14.7229+16.8335i) * in .* absIn(:,2) + (61.6423-76.9171i) * in .* absIn(:,3) + ...
        (-145.139+184.765i) * in .* absIn(:,4) + (190.61-239.371i)* in .* absIn(:,5) + ...
        (-130.184+158.957i) * in .* absIn(:,6) + (36.0047-42.5192i) * in .* absIn(:,7);

end

function out = paMemorylessGaN2Dot1GHz(in)
    % 2.1 GHz GaN
    absIn = abs(in).^(2*(1:4));
    out = (0.999952-0.00981788i) * in + (-0.0618171+0.118845i) * in .* absIn(:,1) + ...
        (-1.69917-0.464933i) * in .* absIn(:,2) + (3.27962+0.829737i) * in .* absIn(:,3) + ...
        (-1.80821-0.454331i) * in .* absIn(:,4);

end

function out = paMemorylessCMOS28GHz(in)
    % 28 GHz CMOS
    absIn = abs(in).^(2*(1:7));
    out = (0.491576+0.870835i) * in + (-1.26213+0.242689i) * in .* absIn(:,1) + ...
        (7.11693+5.14105i) * in .* absIn(:,2) + (-30.7048-53.4924i) * in .* absIn(:,3) + ...
        (73.8814+169.146i) * in .* absIn(:,4) + (-96.7955-253.635i)* in .* absIn(:,5) + ...
        (65.0665+185.434i) * in .* absIn(:,6) + (-17.5838-53.1786i) * in .* absIn(:,7);

end

function out = paMemorylessGaN28GHz(in)
    % 28 GHz GaN
    absIn = abs(in).^(2*(1:5));
    out = (-0.334697-0.942326i) * in + (0.89015-0.72633i) * in .* absIn(:,1) + ...
        (-2.58056+4.81215i) * in .* absIn(:,2) + (4.81548-9.54837i) * in .* absIn(:,3) + ...
        (-4.41452+8.63164i) * in .* absIn(:,4) + (1.54271-2.94034i)* in .* absIn(:,5);

end

function paChar = getDefaultLookup
    % The operating specification for the LDMOS-based Doherty
    % amplifier are:
    % * A frequency of 2110 MHz
    % * A peak power of 300 W
    % * A small signal gain of 61 dB
    % Each row in HAV08_Table specifies Pin (dBm), gain (dB), phase
    % shift (degrees) as derived from figure 4 of Hammi, Oualid, et
    % al. "Power amplifiers' model assessment and memory effects
    % intensity quantification using memoryless post-compensation
    % technique." IEEE Transactions on Microwave Theory and
    % Techniques 56.12 (2008): 3170-3179.

    HAV08_Table =...
        [-35,60.53,0.01;
        -34,60.53,0.01;
        -33,60.53,0.08;
        -32,60.54,0.08;
        -31,60.55,0.1;
        -30,60.56,0.08;
        -29,60.57,0.14;
        -28,60.59,0.19;
        -27,60.6,0.23;
        -26,60.64,0.21;
        -25,60.69,0.28;
        -24,60.76,0.21;
        -23,60.85,0.12;
        -22,60.97,0.08;
        -21,61.12,-0.13;
        -20,61.31,-0.44;
        -19,61.52,-0.94;
        -18,61.76,-1.59;
        -17,62.01,-2.73;
        -16,62.25,-4.31;
        -15,62.47,-6.85;
        -14,62.56,-9.82;
        -13,62.47,-12.29;
        -12,62.31,-13.82;
        -11,62.2,-15.03;
        -10,62.15,-16.27;
        -9,62,-18.05;
        -8,61.53,-20.21;
        -7,60.93,-23.38;
        -6,60.2,-26.64;
        -5,59.38,-28.75];
    % Convert the second column of the HAV08_Table from gain to
    % Pout for use by the memoryless nonlinearity System object.
    paChar = HAV08_Table;
    paChar(:,2) = paChar(:,1) + paChar(:,2);
end

function out = getDefaultCoefficients
    % The 2.44 GHz memory polynomial model defined in TR 38.803
    % Appendix A. Memory-polynomial depth is 5 and
    % memory-polynomial degree is 5. Rows in the output corresponds
    % to memory depth.
    out = [20.0875+0.4240i -6.3792-0.5507i 0.5809+0.0644i 1.6619+0.1040i -0.3561-0.1033i; ...
        -59.8327-34.7815i -2.4805+0.9344i 4.2741+0.7696i -2.0014-2.3785i -1.2566+1.0495i; ...
        3.2738e2+8.4121e2i 4.4019e2-3.0714e1i -3.5935e2-9.9152e0i 1.6961e2+7.3829e1i -4.1661-21.1090i; ...
        -1.6352e3-5.5757e3i -2.5782e3+3.3332e2i 1.9915e3-1.4479e2i -9.0167e2-5.4617e2i -93.1907+14.2774i; ...
        2.3022e3+1.2348e4i 4.6476e3-1.4477e3i -2.9998e3+1.6071e3i 9.1856e2+9.8066e2i 8.2544e2+6.1424e2i].';
end