NR NTN PDSCH Throughput

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)

displaySimulationInformation = true;

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

scaleFactor = []; % Amplitude scaling, in dB

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;

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;

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

waveformInfo = nrOFDMInfo(simParameters.Carrier);

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

% 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

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;

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