NR NTN PDSCH Throughput
This example shows how to measure the physical downlink shared channel (PDSCH) throughput of a 5G New Radio (NR) link in a non-terrestrial network (NTN) channel, as defined by the 3GPP NR standard. The example implements the PDSCH and downlink shared channel (DL-SCH). The transmitter model includes PDSCH demodulation reference signals (DM-RS) and PDSCH phase tracking reference signals (PT-RS). The example supports NTN narrowband and NTN tapped delay line (TDL) propagation channels.
Introduction
This example measures the PDSCH throughput of a 5G link, as defined by the 3GPP NR standard [1], [2], [3], [4].
The example models these 5G NR features:
DL-SCH transport channel coding
Multiple codewords, dependent on the number of layers
PDSCH, PDSCH DM-RS, and PDSCH PT-RS generation
Variable subcarrier spacing and frame numerologies
Normal and extended cyclic prefix
NTN narrowband and NTN TDL propagation channel models
Other features of the simulation are:
PDSCH precoding using singular value decomposition (SVD)
Cyclic prefix orthogonal frequency division multiplexing (CP-OFDM) modulation
Slot-wise and non-slot-wise PDSCH and DM-RS mapping
Timing synchronization and channel estimation
A single bandwidth part (BWP) across the whole carrier
The figure shows the implemented processing chain. For clarity, DM-RS and PT-RS generation are omitted.
For a more detailed explanation of the steps implemented in this example, see Model 5G NR Communication Links (5G Toolbox) and DL-SCH and PDSCH Transmit and Receive Processing Chain (5G Toolbox).
This example supports wideband and subband precoding. You determine the precoding matrix by using SVD and averaging the channel estimate across all PDSCH PRBs either in the allocation (for wideband precoding) or in the subband.
To reduce the total simulation time, you can use Parallel Computing Toolbox™ to execute the SNR points of the SNR loop in parallel.
Simulation Length and SNR Points
Set the length of the simulation in terms of the number of 10 ms frames. By default, the example uses 2 frames, but, a large number of 10 ms frames is necessary to produce meaningful throughput results. Set the SNR points to simulate. The SNR for each layer is defined per resource element (RE), and it includes the effect of signal and noise across all antennas. For an explanation of the SNR definition that this example uses, see SNR Definition Used in Link Simulations (5G Toolbox).
simParameters = struct; % Clear simParameters variable to contain % all key simulation parameters simParameters.NFrames = 2; % Number of 10 ms frames simParameters.SNRIn = 5:2:15; % SNR range (dB)
Set the displaySimulationInformation
variable to true
to display information about the throughput simulation at each SNR point.
displaySimulationInformation =
true;
Doppler Compensation Configuration
The example supports two Doppler compensation configurations: one at transmitter end and the other at receiver end. For compensation at transmitter end, use the DopplerPreCompensator
field. Set the DopplerPreCompensator
field to true to account for Doppler due to satellite movement by applying Doppler pre-compensation to the transmitted waveform. For compensation at receiver end, use the RxDopplerCompensator
field. Set the RxDopplerCompensator
field to true to estimate and compensate the Doppler shift for the received waveform.
simParameters.DopplerPreCompensator =true; simParameters.RxDopplerCompensator =
false;
Carrier and PDSCH Configuration
Set the key parameters of the simulation. These parameters include:
Bandwidth in resource blocks (RBs)
Subcarrier spacing (SCS) in kHz: 15, 30, 60, or 120
Cyclic prefix length (CP): normal or extended
Cell identity
Number of transmit and receive antennas
Also create a substructure containing the DL-SCH and PDSCH parameters, including:
Target code rate
Allocated resource blocks (PRBSet)
Modulation scheme: 'QPSK', '16QAM', '64QAM', or '256QAM'
Number of layers
PDSCH mapping type
DM-RS configuration parameters
PT-RS configuration parameters
% Set waveform type and PDSCH numerology (SCS and CP type) simParameters.Carrier = nrCarrierConfig; simParameters.Carrier.SubcarrierSpacing =30; simParameters.Carrier.CyclicPrefix =
'Normal'; % Bandwidth in number of RBs (11 RBs at 30 kHz SCS for 5 MHz bandwidth) simParameters.Carrier.NSizeGrid = 11; % Physical layer cell identity simParameters.Carrier.NCellID =
1; % PDSCH/DL-SCH parameters % This PDSCH definition is the basis for all PDSCH transmissions in the % throughput simulation simParameters.PDSCH = nrPDSCHConfig; % This structure is to hold additional simulation parameters for the DL-SCH % and PDSCH simParameters.PDSCHExtension = struct(); % Define PDSCH time-frequency resource allocation per slot to be full grid % (single full grid BWP) % PDSCH PRB allocation simParameters.PDSCH.PRBSet = 0:simParameters.Carrier.NSizeGrid-1; % Starting symbol and number of symbols of each PDSCH allocation simParameters.PDSCH.SymbolAllocation = [0,simParameters.Carrier.SymbolsPerSlot]; simParameters.PDSCH.MappingType =
'A'; % Scrambling identifiers simParameters.PDSCH.NID = simParameters.Carrier.NCellID; simParameters.PDSCH.RNTI =
1; % PDSCH resource block mapping (TS 38.211 Section 7.3.1.6) simParameters.PDSCH.VRBToPRBInterleaving =
0; simParameters.PDSCH.VRBBundleSize =
4; % Define the number of transmission layers to be used simParameters.PDSCH.NumLayers =
1; % Define codeword modulation and target coding rate % The number of codewords is directly dependent on the number of layers so % ensure that layers are set first before getting the codeword number if simParameters.PDSCH.NumCodewords > 1 % Multicodeword transmission (when number of layers being > 4) simParameters.PDSCH.Modulation = {
'16QAM',
'16QAM'}; % Code rate used to calculate transport block sizes simParameters.PDSCHExtension.TargetCodeRate = [490 490]/1024; else simParameters.PDSCH.Modulation =
'16QAM'; % Code rate used to calculate transport block size simParameters.PDSCHExtension.TargetCodeRate = 490/1024; end % DM-RS and antenna port configuration (TS 38.211 Section 7.4.1.1) simParameters.PDSCH.DMRS.DMRSPortSet = 0:simParameters.PDSCH.NumLayers-1; simParameters.PDSCH.DMRS.DMRSTypeAPosition =
2; simParameters.PDSCH.DMRS.DMRSLength =
1; simParameters.PDSCH.DMRS.DMRSAdditionalPosition =
2; simParameters.PDSCH.DMRS.DMRSConfigurationType =
2; simParameters.PDSCH.DMRS.NumCDMGroupsWithoutData =
1; simParameters.PDSCH.DMRS.NIDNSCID =
1; simParameters.PDSCH.DMRS.NSCID =
0; % PT-RS configuration (TS 38.211 Section 7.4.1.2) simParameters.PDSCH.EnablePTRS =
0; simParameters.PDSCH.PTRS.TimeDensity =
1; simParameters.PDSCH.PTRS.FrequencyDensity =
2; simParameters.PDSCH.PTRS.REOffset =
'00'; % PT-RS antenna port, subset of DM-RS port set. Empty corresponds to lowest % DM-RS port number simParameters.PDSCH.PTRS.PTRSPortSet = []; % Reserved PRB patterns, if required (for CORESETs, forward compatibility etc) simParameters.PDSCH.ReservedPRB{1}.SymbolSet = []; % Reserved PDSCH symbols simParameters.PDSCH.ReservedPRB{1}.PRBSet = []; % Reserved PDSCH PRBs simParameters.PDSCH.ReservedPRB{1}.Period = []; % Periodicity of reserved resources % Additional simulation and DL-SCH related parameters % PDSCH PRB bundling (TS 38.214 Section 5.1.2.3) simParameters.PDSCHExtension.PRGBundleSize = []; % 2, 4, or [] to signify "wideband" % Rate matching or transport block size (TBS) parameters % Set PDSCH rate matching overhead for TBS (Xoh) to 6 when PT-RS is enabled, otherwise 0 simParameters.PDSCHExtension.XOverhead = 6*simParameters.PDSCH.EnablePTRS; % Hybrid automatic repeat request (HARQ) parameters % Number of parallel HARQ processes to use simParameters.PDSCHExtension.NHARQProcesses = 1; % Enable retransmissions for each process, using redundancy version (RV) sequence [0,2,3,1] simParameters.PDSCHExtension.EnableHARQ = false; % LDPC decoder parameters % Available algorithms: 'Belief propagation', 'Layered belief propagation', % 'Normalized min-sum', 'Offset min-sum' simParameters.PDSCHExtension.LDPCDecodingAlgorithm =
'Normalized min-sum'; simParameters.PDSCHExtension.MaximumLDPCIterationCount = 6; % Define the overall transmission antenna geometry at end-points % For NTN narrowband channel, only single-input-single-output (SISO) % transmission is allowed % Number of PDSCH transmission antennas (1,2,4,8,16,32,64,128,256,512,1024) >= NumLayers simParameters.NumTransmitAntennas = 1; if simParameters.PDSCH.NumCodewords > 1 % Multi-codeword transmission % Number of UE receive antennas (even number >= NumLayers) simParameters.NumReceiveAntennas = 8; else % Number of UE receive antennas (1 or even number >= NumLayers) simParameters.NumReceiveAntennas = 1; end
Get information about the baseband waveform after the OFDM modulation step.
waveformInfo = nrOFDMInfo(simParameters.Carrier);
Propagation Channel Model Construction
Create the channel model object for the simulation. Both the NTN narrowband and NTN TDL channel models are supported [5], [6]. For more information on how to model NTN narrowband and NTN TDL channels, see Model NR NTN Channel.
% Define the general NTN propagation channel parameters % Set the NTN channel type to Narrowband for an NTN narrowband channel and % set the NTN channel type to TDL for a NTN TDL channel. simParameters.NTNChannelType ='Narrowband'; % Set the parameters common to both NTN narrowband and NTN TDL channels simParameters.CarrierFrequency = 2e9; % Carrier frequency (in Hz) simParameters.ElevationAngle = 50; % Elevation angle (in degrees) simParameters.MobileSpeed = 3*1000/3600; % Speed of mobile terminal (in m/s) simParameters.SatelliteAltitude = 600000; % Satellite altitude (in m) simParameters.SatelliteSpeed = 7562.2; % Satellite speed (in m/s) simParameters.SampleRate = waveformInfo.SampleRate; simParameters.RandomStream = 'mt19937ar with seed'; simParameters.Seed = 73; % Set the following fields for NTN narrowband channel if strcmpi(simParameters.NTNChannelType,'Narrowband') simParameters.Environment = 'Urban'; simParameters.AzimuthOrientation = 0; end % Set the following fields for NTN TDL channel if strcmpi(simParameters.NTNChannelType,'TDL') simParameters.DelayProfile = 'NTN-TDL-A'; simParameters.DelaySpread = 30e-9; end % Cross-check the PDSCH layering against the channel geometry validateNumLayers(simParameters); % Set up the NTN channel ntnChan = HelperSetupNTNChannel(simParameters); % Get the maximum number of delayed samples due to a channel multipath % component. The maximum number of delayed samples is calculated from the % channel path with the maximum delay and the implementation delay of the % channel filter. This number of delay samples is required later to flush % the channel filter to obtain the received signal. chInfo = info(ntnChan.BaseChannel); maxChDelay = ceil(max(chInfo.PathDelays*ntnChan.BaseChannel.SampleRate)) + ... chInfo.ChannelFilterDelay;
Processing Loop
To determine the throughput at each SNR point, analyze the PDSCH data for each transmission instance using these steps.
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.
Generate the resource grid — The
nrDLSCH
(5G Toolbox) System object performs transport channel coding. The object operates on the input transport block. ThenrPDSCH
(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.Generate the waveform — The
nrOFDMModulate
(5G Toolbox) function OFDM-modulates the generated grid to get the time-domain waveform.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.
Model and apply a noisy channel — Pass the generated waveform through an NTN narrowband or NTN TDL fading channel to get the faded waveform. Add additive white Gaussian noise (AWGN) to the faded waveform. Define the SNR for each layer based on the RE and receive antenna. For an explanation of the SNR definition that this example uses, see SNR Definition Used in Link Simulations (5G Toolbox).
Apply Doppler compensation — Estimate the Doppler shift in the received waveform and compensate the Doppler shift.
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.Perform channel estimation — For channel estimation, PDSCH DM-RS is used.
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.Calculate precoding matrix — Use SVD to generate the precoding matrix W for the next transmission.
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.Decode the DL-SCH — Pass the decoded soft bits through the
nrDLSCHDecoder
(5G Toolbox) System object. The object decodes the codeword and returns the block cyclic redundancy check (CRC) error. Update the HARQ process with the CRC error. This example determines the throughput of the PDSCH link using the CRC error.
% Array to store the maximum throughput for all SNR points maxThroughput = zeros(length(simParameters.SNRIn),1); % Array to store the simulation throughput for all SNR points simThroughput = zeros(length(simParameters.SNRIn),1); % Set up RV sequence for all HARQ processes if simParameters.PDSCHExtension.EnableHARQ % In the final report of RAN WG1 meeting #91 (R1-1719301), it was % observed in R1-1717405 that if performance is the priority, [0 2 3 1] % should be used. If self-decodability is the priority, it should be % taken into account that the upper limit of the code rate at which % each RV is self-decodable is in the following order: 0>3>2>1 rvSeq = [0 2 3 1]; else % In case of HARQ disabled, RV is set to 0 rvSeq = 0; end % Create DL-SCH encoder System object to perform transport channel encoding encodeDLSCH = nrDLSCH; encodeDLSCH.MultipleHARQProcesses = true; encodeDLSCH.TargetCodeRate = simParameters.PDSCHExtension.TargetCodeRate; % Create DL-SCH decoder System object to perform transport channel decoding decodeDLSCH = nrDLSCHDecoder; decodeDLSCH.MultipleHARQProcesses = true; decodeDLSCH.TargetCodeRate = simParameters.PDSCHExtension.TargetCodeRate; decodeDLSCH.LDPCDecodingAlgorithm = simParameters.PDSCHExtension.LDPCDecodingAlgorithm; decodeDLSCH.MaximumLDPCIterationCount = ... simParameters.PDSCHExtension.MaximumLDPCIterationCount; % Set a threshold value to detect the valid OFDM symbol boundary. For a % SISO case, a threshold of 0.48 can be used to have probability of % incorrect boundary detection around 0.01. Use 0 to avoid thresholding % logic. dtxThresold = 0.48; % Use an offset to account for the common delay. The example, by default, % does not introduce any common delay and only passes through the channel. sampleDelayOffset = 0; % Number of samples % Set this flag to true, to use the shift value estimated in first slot % directly for the consecutive slots. When set to false, the shift is % calculated for each slot, considering the range of shift values to be % whole cyclic prefix length. This is used in the estimation of integer % Doppler shift. usePreviousShift = false; % Processing loop for snrIdx = 1:numel(simParameters.SNRIn) % Comment out for parallel computing % parfor snrIdx = 1:numel(simParameters.SNRIn) % Uncomment for parallel computing % To reduce the total simulation time, you can execute this loop in % parallel by using Parallel Computing Toolbox features. Comment % out the for-loop statement and uncomment the parfor-loop statement. % If Parallel Computing Toolbox is not installed, parfor-loop defaults % to a for-loop statement. Because the parfor-loop iterations are % executed in parallel in a nondeterministic order, the simulation % information displayed for each SNR point can be intertwined. To % switch off the simulation information display, set the % displaySimulationInformation variable (defined earlier in this % example) to false. % Reset the random number generator so that each SNR point experiences % the same noise realization rng('default'); % Make copies of the simulation-level parameter structures so that they % are not Parallel Computing Toolbox broadcast variables when using parfor simLocal = simParameters; waveinfoLocal = waveformInfo; ntnChanLocal = ntnChan; % Make copies of channel-level parameters to simplify subsequent % parameter referencing carrier = simLocal.Carrier; pdsch = simLocal.PDSCH; pdschextra = simLocal.PDSCHExtension; % Copy of the decoder handle to help Parallel Computing Toolbox % classification decodeDLSCHLocal = decodeDLSCH; decodeDLSCHLocal.reset(); % Reset decoder at the start of each SNR point pathFilters = []; thres = dtxThresold; sampleOffset = sampleDelayOffset; usePrevShift = usePreviousShift; % SNR value SNRdB = simLocal.SNRIn(snrIdx); % Specify the order in which we cycle through the HARQ process % identifiers harqSequence = 0:pdschextra.NHARQProcesses-1; % Initialize the state of all HARQ processes harqEntity = HARQEntity(harqSequence,rvSeq,pdsch.NumCodewords); % Reset the channel so that each SNR point experiences the same channel % realization reset(ntnChanLocal.BaseChannel); reset(ntnChanLocal.ChannelFilter); % Total number of slots in the simulation period NSlots = simLocal.NFrames*carrier.SlotsPerFrame; % Obtain a precoding matrix (wtx) to use in the transmission of the % first transport block [estChannelGrid,sampleTimes] = getInitialChannelEstimate(... carrier,simLocal.NumTransmitAntennas,ntnChanLocal); newWtx = getPrecodingMatrix(carrier,pdsch,estChannelGrid,pdschextra.PRGBundleSize); % Timing offset, updated in every slot when the correlation is strong offset = 0; tDoppler = sampleTimes(end); shiftOut = 0; % Loop over the entire waveform length for nslot = 0:NSlots-1 % Update carrier slot number to account for new slot transmission carrier.NSlot = nslot; % Calculate the transport block sizes for the transmission in the slot [pdschIndices,pdschIndicesInfo] = nrPDSCHIndices(carrier,pdsch); trBlkSizes = nrTBS(pdsch.Modulation,pdsch.NumLayers,numel(pdsch.PRBSet),... pdschIndicesInfo.NREPerPRB,pdschextra.TargetCodeRate,pdschextra.XOverhead); % Set transport block depending on the HARQ process for cwIdx = 1:pdsch.NumCodewords % Create a new DL-SCH transport block for new data in the % current process if harqEntity.NewData(cwIdx) trBlk = randi([0 1],trBlkSizes(cwIdx),1); setTransportBlock(encodeDLSCH,trBlk,cwIdx-1,harqEntity.HARQProcessID); % Flush decoder soft buffer explicitly for any new data % because of previous RV sequence time out if harqEntity.SequenceTimeout(cwIdx) resetSoftBuffer(decodeDLSCHLocal,cwIdx-1,harqEntity.HARQProcessID); end end end % Encode the DL-SCH transport blocks codedTrBlocks = encodeDLSCH(pdsch.Modulation,pdsch.NumLayers, ... pdschIndicesInfo.G,harqEntity.RedundancyVersion,harqEntity.HARQProcessID); % Get precoding matrix (wtx) calculated in previous slot wtx = newWtx; % Perform PDSCH modulation pdschSymbols = nrPDSCH(carrier,pdsch,codedTrBlocks); % Create resource grid associated with PDSCH transmission antennas pdschGrid = nrResourceGrid(carrier,simLocal.NumTransmitAntennas); % Perform implementation-specific PDSCH MIMO precoding and mapping [pdschAntSymbols,pdschAntIndices] = hPRGPrecode(size(pdschGrid), ... carrier.NStartGrid,pdschSymbols,pdschIndices,wtx); pdschGrid(pdschAntIndices) = pdschAntSymbols; % Perform implementation-specific PDSCH DM-RS MIMO precoding and % mapping dmrsSymbols = nrPDSCHDMRS(carrier,pdsch); dmrsIndices = nrPDSCHDMRSIndices(carrier,pdsch); [dmrsAntSymbols,dmrsAntIndices] = hPRGPrecode(size(pdschGrid), ... carrier.NStartGrid,dmrsSymbols,dmrsIndices,wtx); pdschGrid(dmrsAntIndices) = dmrsAntSymbols; % Perform implementation-specific PDSCH PT-RS MIMO precoding and % mapping ptrsSymbols = nrPDSCHPTRS(carrier,pdsch); ptrsIndices = nrPDSCHPTRSIndices(carrier,pdsch); [ptrsAntSymbols,ptrsAntIndices] = hPRGPrecode(size(pdschGrid), ... carrier.NStartGrid,ptrsSymbols,ptrsIndices,wtx); pdschGrid(ptrsAntIndices) = ptrsAntSymbols; % Perform OFDM modulation txWaveform = nrOFDMModulate(carrier,pdschGrid); % Pass data through the channel model. Append zeros at the end of % the transmitted waveform to flush the channel content. These % zeros take into account any delay introduced in the channel. This % delay is a combination of the multipath delay and implementation % delay. This value can change depending on the sampling rate, % delay profile, and delay spread. Also apply Doppler % pre-compensation. txWaveform = [txWaveform; zeros(maxChDelay,size(txWaveform,2))]; %#ok<AGROW> [txWaveform,tDoppler] = compensateDopplerShift(... txWaveform,ntnChanLocal.BaseChannel.SampleRate, ... ntnChanLocal.SatelliteDopplerShift, ... simLocal.DopplerPreCompensator,tDoppler); [rxWaveform,pathGains] = HelperGenerateNTNChannel(... ntnChanLocal,txWaveform); % Add AWGN to the received time-domain waveform. Normalize the % noise power by the size of the inverse fast Fourier transform % (IFFT) used in OFDM modulation, because the OFDM modulator % applies this normalization to the transmitted waveform. Also, % normalize the noise power by the number of receive antennas, % because the default behavior of the channel model is to apply % this normalization to the received waveform. SNR = 10^(SNRdB/10); noiseScaling = ... 1/sqrt(2.0*simLocal.NumReceiveAntennas*double(waveinfoLocal.Nfft)*SNR); noise = noiseScaling*complex(randn(size(rxWaveform)),randn(size(rxWaveform))); rxWaveform = rxWaveform + noise; % Perform fractional Doppler frequency shift estimation and % compensation. Use the cyclic prefix in the OFDM waveform to % compute the fractional Doppler shift. [fractionalDopplerShift,detFlag] = estimateFractionalDopplerShift( ... rxWaveform,carrier.SubcarrierSpacing,waveinfoLocal.Nfft, ... waveinfoLocal.CyclicPrefixLengths(2),thres, ... simLocal.RxDopplerCompensator); rxWaveform = compensateDopplerShift(rxWaveform,waveinfoLocal.SampleRate, ... fractionalDopplerShift,simLocal.RxDopplerCompensator); % Perform integer Doppler frequency shift estimation and % compensation. Use the demodulation reference signals to compute % the integer Doppler shift. [integerDopplerShift,shiftOut] = estimateIntegerDopplerShift(carrier,rxWaveform, ... dmrsIndices,dmrsSymbols,sampleOffset,usePrevShift,shiftOut-sampleOffset, ... (simLocal.RxDopplerCompensator && detFlag)); rxWaveform = compensateDopplerShift(rxWaveform,waveinfoLocal.SampleRate, ... integerDopplerShift,simLocal.RxDopplerCompensator); % For timing synchronization, correlate the received waveform with % the PDSCH DM-RS to give timing offset estimate t and correlation % magnitude mag. The function hSkipWeakTimingOffset is used to % update the receiver timing offset. If the correlation peak in mag % is weak, the current timing estimate t is ignored and the % previous estimate offset is used. [t,mag] = nrTimingEstimate(carrier,rxWaveform,dmrsIndices,dmrsSymbols); offset = hSkipWeakTimingOffset(offset,t,mag); % Display a warning if the estimated timing offset exceeds the % maximum channel delay if offset > maxChDelay warning(['Estimated timing offset (%d) is greater than the maximum' ... ' channel delay (%d). This might result in a decoding failure.' ... ' The estimated timing offset might be greater than the maximum' ... ' channel delay because of low SNR or not enough DM-RS symbols to' ... ' synchronize successfully.'],offset,maxChDelay); end rxWaveform = rxWaveform(1+offset:end,:); % Perform OFDM demodulation on the received data to recreate the % resource grid. Include zero padding in the event that practical % synchronization results in an incomplete slot being demodulated. rxGrid = nrOFDMDemodulate(carrier,rxWaveform); [K,L,R] = size(rxGrid); if (L < carrier.SymbolsPerSlot) rxGrid = cat(2,rxGrid,zeros(K,carrier.SymbolsPerSlot-L,R)); end % Perform least squares channel estimation between the received % grid and each transmission layer, using the PDSCH DM-RS for each % layer. This channel estimate includes the effect of transmitter % precoding. [estChannelGrid,noiseEst] = nrChannelEstimate(carrier,rxGrid,... dmrsIndices,dmrsSymbols,'CDMLengths',pdsch.DMRS.CDMLengths); % Get PDSCH REs from the received grid and estimated channel grid [pdschRx,pdschHest] = nrExtractResources(... pdschIndices,rxGrid,estChannelGrid); % Remove precoding from estChannelGrid prior to precoding % matrix calculation estChannelGridPorts = precodeChannelEstimate(... carrier,estChannelGrid,conj(wtx)); % Get precoding matrix for next slot newWtx = getPrecodingMatrix(... carrier,pdsch,estChannelGridPorts,pdschextra.PRGBundleSize); % Perform equalization [pdschEq,csi] = nrEqualizeMMSE(pdschRx,pdschHest,noiseEst); % Common phase error (CPE) compensation if ~isempty(ptrsIndices) % Initialize temporary grid to store equalized symbols tempGrid = nrResourceGrid(carrier,pdsch.NumLayers); % Extract PT-RS symbols from received grid and estimated % channel grid [ptrsRx,ptrsHest,~,~,~,ptrsLayerIndices] = ... nrExtractResources(ptrsIndices,rxGrid,estChannelGrid,tempGrid); % Equalize PT-RS symbols and map them to tempGrid ptrsEq = nrEqualizeMMSE(ptrsRx,ptrsHest,noiseEst); tempGrid(ptrsLayerIndices) = ptrsEq; % Estimate the residual channel at the PT-RS locations in % tempGrid cpe = nrChannelEstimate(tempGrid,ptrsIndices,ptrsSymbols); % Sum estimates across subcarriers, receive antennas, and % layers. Then, get the CPE by taking the angle of the % resultant sum cpe = angle(sum(cpe,[1 3 4])); % Map the equalized PDSCH symbols to tempGrid tempGrid(pdschIndices) = pdschEq; % Correct CPE in each OFDM symbol within the range of reference % PT-RS OFDM symbols symLoc = ... pdschIndicesInfo.PTRSSymbolSet(1)+1:pdschIndicesInfo.PTRSSymbolSet(end)+1; tempGrid(:,symLoc,:) = tempGrid(:,symLoc,:).*exp(-1i*cpe(symLoc)); % Extract PDSCH symbols pdschEq = tempGrid(pdschIndices); end % Decode PDSCH symbols [dlschLLRs,rxSymbols] = nrPDSCHDecode(carrier,pdsch,pdschEq,noiseEst); % Scale the decoded log-likelihood ratios (LLRs) by channel state % information (CSI) csi = nrLayerDemap(csi); % CSI layer demapping for cwIdx = 1:pdsch.NumCodewords Qm = length(dlschLLRs{cwIdx})/length(rxSymbols{cwIdx}); % Bits per symbol csi{cwIdx} = repmat(csi{cwIdx}.',Qm,1); % Expand by each bit % per symbol dlschLLRs{cwIdx} = dlschLLRs{cwIdx} .* csi{cwIdx}(:); % Scale by CSI end % Decode the DL-SCH transport channel decodeDLSCHLocal.TransportBlockLength = trBlkSizes; [decbits,blkerr] = decodeDLSCHLocal(dlschLLRs,pdsch.Modulation,... pdsch.NumLayers,harqEntity.RedundancyVersion,harqEntity.HARQProcessID); % Store values to calculate throughput simThroughput(snrIdx) = simThroughput(snrIdx) + sum(~blkerr .* trBlkSizes); maxThroughput(snrIdx) = maxThroughput(snrIdx) + sum(trBlkSizes); % Update current process with CRC error and advance to next process updateAndAdvance(harqEntity,blkerr,trBlkSizes,pdschIndicesInfo.G); end % Display the results if displaySimulationInformation == 1 fprintf('\nThroughput(Mbps) for %d frame(s) at SNR %d dB: %.4f\n',... simLocal.NFrames,SNRdB,1e-6*simThroughput(snrIdx)/(simLocal.NFrames*10e-3)); fprintf('Throughput(%%) for %d frame(s) at SNR %d dB: %.4f\n',... simLocal.NFrames,SNRdB,simThroughput(snrIdx)*100/maxThroughput(snrIdx)); end end
Throughput(Mbps) for 2 frame(s) at SNR 5 dB: 0.0000
Throughput(%) for 2 frame(s) at SNR 5 dB: 0.0000
Throughput(Mbps) for 2 frame(s) at SNR 7 dB: 0.0000
Throughput(%) for 2 frame(s) at SNR 7 dB: 0.0000
Throughput(Mbps) for 2 frame(s) at SNR 9 dB: 6.2308
Throughput(%) for 2 frame(s) at SNR 9 dB: 92.5000
Throughput(Mbps) for 2 frame(s) at SNR 11 dB: 6.7360
Throughput(%) for 2 frame(s) at SNR 11 dB: 100.0000
Throughput(Mbps) for 2 frame(s) at SNR 13 dB: 6.7360
Throughput(%) for 2 frame(s) at SNR 13 dB: 100.0000
Throughput(Mbps) for 2 frame(s) at SNR 15 dB: 6.7360
Throughput(%) for 2 frame(s) at SNR 15 dB: 100.0000
Results
Display the measured throughput, which is the percentage of the maximum possible throughput of the link given the available resources for data transmission.
figure; plot(simParameters.SNRIn,simThroughput*100./maxThroughput,'o-.') xlabel('SNR (dB)'); ylabel('Throughput (%)'); grid on; title(sprintf('NTN %s (%dx%d) / NRB=%d / SCS=%dkHz', ... simParameters.NTNChannelType,simParameters.NumTransmitAntennas, ... simParameters.NumReceiveAntennas, simParameters.Carrier.NSizeGrid,... simParameters.Carrier.SubcarrierSpacing));
% Bundle key parameters and results into a combined structure for recording
simResults.simParameters = simParameters;
simResults.simThroughput = simThroughput;
simResults.maxThroughput = maxThroughput;
This next figure shows the throughput results obtained by simulating 1000 frames (NFrames = 1000
, SNRIn = 5:2:15
) for a carrier with a 30 kHz SCS occupying a 5 MHz transmission bandwidth. The simulation setup includes the default carrier and PDSCH configuration with an NTN narrowband channel. The line corresponding to the Doppler Pre-compensator is achieved by setting the DopplerPreCompensator
field to true and RxDopplerCompensator
field to false. The line corresponding to the Rx Doppler Compensator is achieved by setting the DopplerPreCompensator
field to false and RxDopplerCompensator
field to true.
Further Exploration
You can use this example to further explore these options:
To analyze the throughput at each SNR point for a different satellite orbit, vary the satellite altitude and satellite speed.
To observe the link performance without any Doppler compensation techniques, set
DopplerPreCompensator
andRxDopplerCompensator
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 andRxDopplerCompensator
field to true.To check the throughput performance of different scenarios, change the carrier numerology and the number of transmit and receive antennas, and set the channel model type to TDL.
To compare the throughput performance in an NTN and terrestrial network, use the
nrTDLChannel
(5G Toolbox) and thenrCDLChannel
(5G Toolbox) channel objects as shown in NR PDSCH Throughput (5G Toolbox).
Appendix
The example uses these helper functions:
HARQEntity — Manage a set of parallel HARQ processes
HelperGenerateNTNChannel — Generate NTN channel
HelperSetupNTNChannel — Set up NTN channel
hPRGPrecode — Precoding for PDSCH precoding resource block group (PRG) bundling
hSkipWeakTimingOffset — Skip timing offset estimates with weak correlation
Selected Bibliography
[1] 3GPP TS 38.211. "NR; Physical channels and modulation." 3rd Generation Partnership Project; Technical Specification Group Radio Access Network.
[2] 3GPP TS 38.212. "NR; Multiplexing and channel coding." 3rd Generation Partnership Project; Technical Specification Group Radio Access Network.
[3] 3GPP TS 38.213. "NR; Physical layer procedures for control." 3rd Generation Partnership Project; Technical Specification Group Radio Access Network.
[4] 3GPP TS 38.214. "NR; Physical layer procedures for data." 3rd Generation Partnership Project; Technical Specification Group Radio Access Network.
[5] 3GPP TR 38.901. "Study on channel model for frequencies from 0.5 to 100 GHz." 3rd Generation Partnership Project; Technical Specification Group Radio Access Network.
[6] 3GPP TR 38.811. "Study on new radio (NR) to support non-terrestrial networks." 3rd Generation Partnership Project; Technical Specification Group Radio Access Network.
[7] 3GPP TR 38.821. "Solutions for NR to support non-terrestrial networks (NTN)." 3rd Generation Partnership Project; Technical Specification Group Radio Access Network.
[8] ITU-R Recommendation P.681-11 (08/2019). "Propagation data required for the design systems in the land mobile-satellite service." P Series; Radio wave propagation.
Local Functions
function validateNumLayers(simParameters) % Validate the number of layers, relative to the antenna geometry numlayers = simParameters.PDSCH.NumLayers; ntxants = simParameters.NumTransmitAntennas; nrxants = simParameters.NumReceiveAntennas; if contains(simParameters.NTNChannelType,'Narrowband','IgnoreCase',true) if (ntxants ~= 1) || (nrxants ~= 1) error(['For NTN narrowband channel, ' ... 'the number of transmit and receive antennas must be 1.']); end end antennaDescription = sprintf(... 'min(NumTransmitAntennas,NumReceiveAntennas) = min(%d,%d) = %d', ... ntxants,nrxants,min(ntxants,nrxants)); if numlayers > min(ntxants,nrxants) error('The number of layers (%d) must satisfy NumLayers <= %s', ... numlayers,antennaDescription); end % Display a warning if the maximum possible rank of the channel equals % the number of layers if (numlayers > 2) && (numlayers == min(ntxants,nrxants)) warning(['The maximum possible rank of the channel, given by %s, is equal to' ... ' NumLayers (%d). This may result in a decoding failure under some channel' ... ' conditions. Try decreasing the number of layers or increasing the channel' ... ' rank (use more transmit or receive antennas).'],antennaDescription, ... numlayers); %#ok<SPWRN> end end function [estChannelGrid,sampleTimes] = getInitialChannelEstimate(... carrier,nTxAnts,ntnChan) % Obtain channel estimate before first transmission. Use this function to % obtain a precoding matrix for the first slot. ofdmInfo = nrOFDMInfo(carrier); chInfo = info(ntnChan.BaseChannel); maxChDelay = ceil(max(chInfo.PathDelays*ntnChan.BaseChannel.SampleRate)) ... + chInfo.ChannelFilterDelay; % Temporary waveform (only needed for the sizes) tmpWaveform = zeros(... (ofdmInfo.SampleRate/1000/carrier.SlotsPerSubframe)+maxChDelay,nTxAnts); % Filter through channel and get the path gains and path filters [~,pathGains,sampleTimes] = HelperGenerateNTNChannel(... ntnChan,tmpWaveform); chanInfo = info(ntnChan.ChannelFilter); pathFilters = chanInfo.ChannelFilterCoefficients.'; % Perfect timing synchronization offset = nrPerfectTimingEstimate(pathGains,pathFilters); % Perfect channel estimate estChannelGrid = nrPerfectChannelEstimate(... carrier,pathGains,pathFilters,offset,sampleTimes); end function wtx = getPrecodingMatrix(carrier,pdsch,hestGrid,prgbundlesize) % Calculate precoding matrices for all precoding resource block groups % (PRGs) in the carrier that overlap with the PDSCH allocation % Maximum common resource block (CRB) addressed by carrier grid maxCRB = carrier.NStartGrid + carrier.NSizeGrid - 1; % PRG size if nargin==4 && ~isempty(prgbundlesize) Pd_BWP = prgbundlesize; else Pd_BWP = maxCRB + 1; end % PRG numbers (1-based) for each RB in the carrier grid NPRG = ceil((maxCRB + 1) / Pd_BWP); prgset = repmat((1:NPRG),Pd_BWP,1); prgset = prgset(carrier.NStartGrid + (1:carrier.NSizeGrid).'); [~,~,R,P] = size(hestGrid); wtx = zeros([pdsch.NumLayers P NPRG]); for i = 1:NPRG % Subcarrier indices within current PRG and within the PDSCH % allocation thisPRG = find(prgset==i) - 1; thisPRG = intersect(thisPRG,pdsch.PRBSet(:) + carrier.NStartGrid,'rows'); prgSc = (1:12)' + 12*thisPRG'; prgSc = prgSc(:); if (~isempty(prgSc)) % Average channel estimate in PRG estAllocGrid = hestGrid(prgSc,:,:,:); Hest = permute(mean(reshape(estAllocGrid,[],R,P)),[2 3 1]); % SVD decomposition [~,~,V] = svd(Hest); wtx(:,:,i) = V(:,1:pdsch.NumLayers).'; end end wtx = wtx / sqrt(pdsch.NumLayers); % Normalize by NumLayers end function estChannelGrid = precodeChannelEstimate(carrier,estChannelGrid,W) % Apply precoding matrix W to the last dimension of the channel estimate [K,L,R,P] = size(estChannelGrid); estChannelGrid = reshape(estChannelGrid,[K*L R P]); estChannelGrid = hPRGPrecode([K L R P],carrier.NStartGrid,estChannelGrid, ... reshape(1:numel(estChannelGrid),[K*L R P]),W); estChannelGrid = reshape(estChannelGrid,K,L,R,[]); end function [loc,wMovSum,pho,bestAnt] = detectOFDMSymbolBoundary(rxWave,nFFT,cpLen,thres) % Detect OFDM symbol boundary by calculating correlation of cyclic prefix % Capture the dimensions of received waveform [NSamples,R] = size(rxWave); % Append zeros of length nFFT across each receive antenna to the % received waveform waveformZeroPadded = [rxWave;zeros(nFFT,R)]; % Get the portion of waveform till the last nFFT samples wavePortion1 = waveformZeroPadded(1:end-nFFT,:); % Get the portion of waveform delayed by nFFT wavePortion2 = waveformZeroPadded(1+nFFT:end,:); % Get the energy of each sample in both the waveform portions eWavePortion1 = abs(wavePortion1).^2; eWavePortion2 = abs(wavePortion2).^2; % Initialize the temporary variables wMovSum = zeros([NSamples R]); wEnergyPortion1 = zeros([NSamples+cpLen-1 R]); wEnergyPortion2 = wEnergyPortion1; % Perform correlation for each sample with the sample delayed by nFFT waveCorr = wavePortion1.*conj(wavePortion2); % Calculate the moving sum value for each cpLen samples, across each % receive antenna oneVec = ones(cpLen,1); for i = 1:R wConv = conv(waveCorr(:,i),oneVec); wMovSum(:,i) = wConv(cpLen:end); wEnergyPortion1(:,i) = conv(eWavePortion1(:,i),oneVec); wEnergyPortion2(:,i) = conv(eWavePortion2(:,i),oneVec); end % Get the normalized correlation value for the moving sum matrix pho = abs(wMovSum)./ ... (eps+sqrt(wEnergyPortion1(cpLen:end,:).*wEnergyPortion2(cpLen:end,:))); % Detect the peak locations in each receive antenna based on the % threshold. These peak locations correspond to the starting location % of each OFDM symbol in the received waveform. loc = cell(R,1); m = zeros(R,1); phoFactor = ceil(NSamples/nFFT); phoExt = [pho; -1*ones(phoFactor*nFFT - NSamples,R)]; for col = 1:R p1 = reshape(phoExt(:,i),[],phoFactor); [pks,locTemp] = max(p1); locTemp = locTemp + (0:phoFactor-1).*nFFT; indicesToConsider = pks>=thres; loc{col} = locTemp(indicesToConsider); m(col) = max(pks); end bestAnt = find(m == max(m)); end function [out,detFlag] = estimateFractionalDopplerShift(rxWave,scs, ... nFFT,cpLen,thres,flag) % Estimate the fractional Doppler shift using cyclic prefix if flag % Detect the OFDM boundary locations [loc,wMovSum,~,bestAnt] = detectOFDMSymbolBoundary(rxWave, ... nFFT,cpLen,thres); % Get the average correlation value at the peak locations for the % first receive antenna having maximum correlation value wSamples = nan(1,1); if ~isempty(loc{bestAnt(1)}) wSamples(1) = mean(wMovSum(loc{bestAnt(1)},bestAnt(1))); end % Compute the fractional Doppler shift if ~all(isnan(wSamples)) out = -(mean(angle(wSamples),'omitnan')*scs*1e3)/(2*pi); % Flag to indicate that there is at least one OFDM symbol % detected detFlag = 1; else out = 0; detFlag = 0; end else out = 0; detFlag = 0; end end function [out,shiftOut] = estimateIntegerDopplerShift(carrier,rx,refInd, ... refSym,sampleOffset,usePrevShift,shiftIn,flag) % Estimate the integer Doppler shift using demodulation reference signal if flag % Get OFDM information ofdmInfo = nrOFDMInfo(carrier); K = carrier.NSizeGrid*12; % Number of subcarriers L = ofdmInfo.SymbolsPerSlot; % Number of OFDM symbols in slot P = ceil(max(double(refInd(:))/(K*L))); % Number of layers cpLen = ofdmInfo.CyclicPrefixLengths(1); % Highest cyclic prefix length % Range of shift values to be used in integer frequency offset % estimation shiftValues = sampleOffset + shiftIn; if ~(usePrevShift && (shiftIn > 0)) % Update range of shift values such that whole cyclic prefix % length is covered shiftValues = sampleOffset + (1:cpLen); end % Initialize temporary variables shiftLen = length(shiftValues); maxValue = complex(zeros(shiftLen,1)); binIndex = zeros(shiftLen,1); R = size(rx,2); rx1 = [rx; zeros(1*(mod(size(rx,1),2)),R)]; % Append zero, if required rxLen = size(rx1,1); x_wave = zeros([rxLen P]); % Generate reference waveform refGrid = complex(zeros([K L P])); refGrid(refInd) = refSym; refWave = nrOFDMModulate(carrier,refGrid,'Windowing',0); refWave = [refWave; zeros((rxLen-size(refWave,1)),P)]; % Find the fast Fourier transform (FFT) bin corresponding to % maximum correlation value for each shift value for shiftIdx = 1:shiftLen % Use the waveform from the shift value and append zeros tmp = rx1(shiftValues(shiftIdx):end,:); rx = [tmp; zeros(rxLen-size(tmp,1),R)]; % Compute the correlation of received waveform with reference % waveform across different layers and receive antennas for rIdx = 1:R for p = 1:P x_wave(:,rIdx,p) = ... rx(:,rIdx).*conj(refWave(1:length(rx(:,rIdx)),p)); end end % Aggregate the correlated waveform across multiple ports and % compute energy of the resultant for each receive antenna x1 = sum(x_wave,3); x1P = sum(abs(x1).^2); % Find the index of first receive antenna which has maximum % correlation energy idx = find(x1P == max(x1P),1); % Combine the received waveform which have maximum correlation % energy x_wave_combined = sum(x1(:,idx(1)),2); % Compute FFT of the resultant waveform x_fft = fftshift(fft(x_wave_combined)); % Store the value and location of peak [maxValue(shiftIdx),binIndex(shiftIdx)] = max(x_fft); end % FFT bin values fftBinValues = (-rxLen/2:(rxLen/2-1))*(ofdmInfo.SampleRate/rxLen); % Find the shift index that corresponds to the maximum of peak % value of all the shifted waveforms. Use the FFT bin index % corresponding to this maximum shift index. The FFT bin value % corresponding to this bin index is the integer frequency offset. [~,maxId] = max(maxValue); loc = binIndex(maxId); out = fftBinValues(loc); shiftOut = shiftValues(maxId); else out = 0; shiftOut = 1+sampleOffset; end end function [out,t] = compensateDopplerShift(inWave,fs,fdSat,flag,t) % Perform Doppler shift correction t1 = (0:size(inWave,1)-1)'/fs; if nargin == 5 t1 = t1 + t; % Add the sample time offset end if flag out = inWave.*exp(1j*2*pi*(-fdSat)*t1); else out = inWave; end t = t1(end); end