ti_mo
1222 lines
56 KiB
Matlab
Executable File
1222 lines
56 KiB
Matlab
Executable File
classdef HelperNRNTNThroughput
|
||
%HelperNRNTNThroughput Class defining the supporting functions used in
|
||
%the NR NTN PDSCH Throughput example
|
||
%
|
||
% Note: This is an undocumented class and its API and/or
|
||
% functionality may change in subsequent releases.
|
||
|
||
% Copyright 2021-2026 The MathWorks, Inc.
|
||
|
||
methods (Static)
|
||
function validateNumLayers(simParameters)
|
||
% Validate the number of layers, relative to the antenna geometry
|
||
|
||
if isfield(simParameters, 'PDSCH')
|
||
numlayers = simParameters.PDSCH.NumLayers;
|
||
else
|
||
numlayers = simParameters.PUSCH.NumLayers;
|
||
end
|
||
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
|
||
|
||
% For NTN CDL channel, the parameters NumTrasnmitAntennas and
|
||
% NumReceiveAntennas must align with that of TxArraySize and
|
||
% RxArraySize respectively
|
||
if simParameters.NTNChannelType == "CDL"
|
||
numElementsTxArray = prod(simParameters.TxArraySize);
|
||
if numElementsTxArray ~= ntxants
|
||
error(['For NTN CDL channel, the product of all values in TxArraySize (%d) must be equal to ' ...
|
||
'number of transmit antennas (%d)'],numElementsTxArray,ntxants);
|
||
end
|
||
numElementsRxArray = prod(simParameters.RxArraySize);
|
||
if numElementsRxArray ~= nrxants
|
||
error(['For NTN CDL channel, the product of all values in RxArraySize (%d) must be equal to ' ...
|
||
'number of receive antennas (%d)'],numElementsRxArray,nrxants);
|
||
end
|
||
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') || isa(channel,'nrCDLChannel')
|
||
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] = ...
|
||
HelperNRNTNThroughput.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
|
||
|
||
arguments
|
||
carrier
|
||
rx
|
||
refInd
|
||
refSym
|
||
sampleOffset = 0
|
||
usePrevShift = false
|
||
useDiffCorr = true
|
||
shiftIn = 0
|
||
maxOffset = 0
|
||
flag = false
|
||
end
|
||
|
||
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 = HelperNRNTNThroughput.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
|
||
|
||
if flag
|
||
out = frequencyOffset(inWave,fs,-fdSat);
|
||
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,varargin)
|
||
% Perform joint time-frequency synchronization
|
||
% jointTimeFreq(carrier,rx,refInd,refSym,fSearchSpace)
|
||
% jointTimeFreq(carrier,rx,refGrid,fSearchSpace)
|
||
fSearchSpace = varargin{end};
|
||
numFreqVals = length(fSearchSpace);
|
||
peakVal = zeros(numFreqVals,1);
|
||
peakIdx = peakVal;
|
||
ofdmInfo = nrOFDMInfo(carrier);
|
||
fs = ofdmInfo.SampleRate;
|
||
for fIdx = 1:numFreqVals
|
||
rxCorrected = ...
|
||
HelperNRNTNThroughput.compensateDopplerShift( ...
|
||
rx,fs,fSearchSpace(fIdx),true);
|
||
[~,corr] = nrTimingEstimate(carrier,rxCorrected,varargin{1:end-1});
|
||
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
|
||
|
||
function [txWaveform,info] = generatePDSCHWaveform(carrier,pdsch,dlsch,wtx,dt,numSlots)
|
||
arguments
|
||
carrier
|
||
pdsch
|
||
dlsch = struct
|
||
wtx = eye(pdsch.NumLayers)
|
||
dt = "double"
|
||
numSlots = 1
|
||
end
|
||
|
||
% Initialize variables
|
||
nTx = size(wtx,2);
|
||
tmpGrid = nrResourceGrid(carrier,nTx,OutputDataType=dt);
|
||
pdschGrid = repmat(tmpGrid,[1 numSlots 1]);
|
||
refGrid = pdschGrid;
|
||
nSlotSymb = carrier.SymbolsPerSlot;
|
||
initialNSlot = carrier.NSlot;
|
||
numCW = pdsch.NumCodewords;
|
||
|
||
% Process loop for each slot
|
||
for slotIdx = 0:numSlots-1
|
||
[slotGrid,refSlotGrid] = deal(tmpGrid);
|
||
carrier.NSlot = initialNSlot + slotIdx;
|
||
|
||
% Perform PDSCH modulation
|
||
[pdschIndices,pdschIndicesInfo] = nrPDSCHIndices(carrier,pdsch);
|
||
if isempty(fieldnames(dlsch))
|
||
% Update the codeword with valid bits
|
||
cw = cell(1,numCW);
|
||
for i = 1:numCW
|
||
cw{i} = randi([0 1],pdschIndicesInfo.G(i),1);
|
||
end
|
||
else
|
||
% Encode with the inputs provided in the dlsch
|
||
% structure
|
||
cw = dlsch.Encoder(pdsch.Modulation,pdsch.NumLayers, ...
|
||
pdschIndicesInfo.G,dlsch.RedundancyVersion, ...
|
||
dlsch.HARQProcessID);
|
||
end
|
||
if ~isempty(cw)
|
||
pdschSymbols = nrPDSCH(carrier,pdsch,cw,OutputDataType=dt);
|
||
% Perform implementation-specific PDSCH MIMO precoding
|
||
% and mapping
|
||
[pdschAntSymbols,pdschAntIndices] = nrPDSCHPrecode( ...
|
||
carrier,pdschSymbols,pdschIndices,wtx);
|
||
slotGrid(pdschAntIndices) = pdschAntSymbols;
|
||
end
|
||
|
||
% Perform implementation-specific PDSCH DM-RS MIMO
|
||
% precoding and mapping
|
||
dmrsSymbols = nrPDSCHDMRS(carrier,pdsch,OutputDataType=dt);
|
||
dmrsIndices = nrPDSCHDMRSIndices(carrier,pdsch);
|
||
[dmrsAntSymbols,dmrsAntIndices] = nrPDSCHPrecode( ...
|
||
carrier,dmrsSymbols,dmrsIndices,wtx);
|
||
slotGrid(dmrsAntIndices) = dmrsAntSymbols;
|
||
refSlotGrid(dmrsAntIndices) = dmrsAntSymbols;
|
||
|
||
% Perform implementation-specific PDSCH PT-RS MIMO
|
||
% precoding and mapping
|
||
ptrsSymbols = nrPDSCHPTRS(carrier,pdsch,OutputDataType=dt);
|
||
ptrsIndices = nrPDSCHPTRSIndices(carrier,pdsch);
|
||
[ptrsAntSymbols,ptrsAntIndices] = nrPDSCHPrecode( ...
|
||
carrier,ptrsSymbols,ptrsIndices,wtx);
|
||
slotGrid(ptrsAntIndices) = ptrsAntSymbols;
|
||
refSlotGrid(ptrsAntIndices) = ptrsAntSymbols;
|
||
|
||
% Map to the grid containing multiple slots
|
||
symIdx = (nSlotSymb*slotIdx)+1:(nSlotSymb*(slotIdx+1));
|
||
pdschGrid(:,symIdx,:) = slotGrid;
|
||
refGrid(:,symIdx,:) = refSlotGrid;
|
||
end
|
||
|
||
% Perform OFDM modulation
|
||
carrier.NSlot = initialNSlot;
|
||
[txWaveform,info] = nrOFDMModulate(carrier,pdschGrid);
|
||
|
||
% Append the resource grids to output structure
|
||
info.ResourceGrid = pdschGrid;
|
||
info.ReferenceGrid = refGrid;
|
||
end
|
||
|
||
function [slotTimes,symLen] = getSlotTimes(nSymbSlot,sl,fs,nf,dt)
|
||
% Get the slot time
|
||
symLen = cumsum(sl); % Lengths of each OFDM symbol in a slot
|
||
nSubFrames = 10;
|
||
samples = [0 symLen(nSymbSlot:nSymbSlot:end-1)]'./fs;
|
||
subframeTimes = (0:(nf*nSubFrames)-1)*1e-3;
|
||
slotTimes = cast(reshape(samples+subframeTimes,[],1),dt);
|
||
end
|
||
|
||
function info = initializeDelayObjects(simParameters, waveformInfo, varargin)
|
||
%initializeDelayObjects Initialize delay objects for NR NTN links.
|
||
%
|
||
% info = initializeDelayObjects(simParameters, waveformInfo)
|
||
% computes delay and path loss based on simParameters and waveformInfo.
|
||
%
|
||
% info = initializeDelayObjects(simParameters, waveformInfo, delayInSeconds)
|
||
% uses the provided delayInSeconds instead of computing it.
|
||
|
||
% Parse optional input
|
||
if nargin > 2 && ~isempty(varargin{1})
|
||
delayInSeconds = varargin{1};
|
||
computeGeometry = false;
|
||
if simParameters.EnableDelay==0
|
||
delayModel = "None";
|
||
else
|
||
delayModel = simParameters.LinkVariationModel;
|
||
end
|
||
else
|
||
computeGeometry = true;
|
||
delayModel = simParameters.DelayModel;
|
||
end
|
||
|
||
% Compute geometry-based parameters if needed
|
||
if computeGeometry
|
||
c = physconst("lightspeed");
|
||
lambda = c / simParameters.CarrierFrequency;
|
||
slotTimes = HelperNRNTNThroughput.getSlotTimes( ...
|
||
simParameters.Carrier.SymbolsPerSlot, waveformInfo.SymbolLengths, ...
|
||
waveformInfo.SampleRate, simParameters.NFrames, simParameters.DataType);
|
||
if (delayModel == "Time-varying")
|
||
SU = slantRangeCircularOrbit(simParameters.ElevationAngle, ...
|
||
simParameters.SatelliteAltitude, simParameters.MobileAltitude, slotTimes);
|
||
lambda = c / simParameters.CarrierFrequency;
|
||
pathLoss = fspl(double(SU), lambda) .* simParameters.IncludeFreeSpacePathLoss;
|
||
delayInSeconds = SU ./ c;
|
||
else
|
||
SU = slantRangeCircularOrbit(simParameters.ElevationAngle, ...
|
||
simParameters.SatelliteAltitude, simParameters.MobileAltitude);
|
||
pathLoss = fspl(SU, lambda) * simParameters.IncludeFreeSpacePathLoss; % dB
|
||
pathLoss = repmat(pathLoss,1,numel(slotTimes)); % dB
|
||
delayInSeconds = cast(SU ./ c, simParameters.DataType);
|
||
end
|
||
|
||
else
|
||
SU = [];
|
||
pathLoss = [];
|
||
end
|
||
|
||
% Initialize configuration objects for delay modeling
|
||
if delayModel == "None"
|
||
staticDelay = dsp.Delay(Length=0);
|
||
variableFractionalDelay = 0;
|
||
variableIntegerDelay = 0;
|
||
delayInSeconds = 0;
|
||
numVariableFracDelaySamples = 0;
|
||
numVariableIntegSamples = 0;
|
||
maxVarPropDelay = 0;
|
||
elseif delayModel == "Static"
|
||
delayInSamples = delayInSeconds .* waveformInfo.SampleRate;
|
||
integDelaySamples = floor(delayInSamples);
|
||
fracDelaySamples = delayInSamples - integDelaySamples;
|
||
numVariableFracDelaySamples = repmat(fracDelaySamples, 1, ...
|
||
(simParameters.Carrier.SlotsPerFrame * simParameters.NFrames) + 1);
|
||
|
||
staticDelay = dsp.Delay(Length=integDelaySamples);
|
||
variableFractionalDelay = dsp.VariableFractionalDelay( ...
|
||
InterpolationMethod="Farrow", ...
|
||
FarrowSmallDelayAction="Use off-centered kernel", ...
|
||
MaximumDelay=1);
|
||
variableIntegerDelay = 0;
|
||
numVariableIntegSamples = 0;
|
||
maxVarPropDelay = 0;
|
||
elseif delayModel == "Time-varying"
|
||
delayInSamples = delayInSeconds .* waveformInfo.SampleRate;
|
||
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);
|
||
else
|
||
error('Unknown DelayModel: %s', delayModel);
|
||
end
|
||
|
||
% Build output structure
|
||
info = struct;
|
||
info.StaticDelay = staticDelay;
|
||
info.VariableIntegerDelay = variableIntegerDelay;
|
||
info.VariableFractionalDelay = variableFractionalDelay;
|
||
info.MaxVariablePropDelay = maxVarPropDelay;
|
||
info.NumVariableIntegerDelaySamples = numVariableIntegSamples;
|
||
info.NumVariableFractionalDelaySamples = numVariableFracDelaySamples;
|
||
info.DelayInSeconds = delayInSeconds;
|
||
|
||
% Add path loss and slant distance only if computed
|
||
if ~isempty(SU)
|
||
info.PathLoss = pathLoss;
|
||
info.SlantDistance = SU;
|
||
end
|
||
|
||
end
|
||
|
||
function [hpa,hpaDelay,paInputScaleFactor] = initializePA( ...
|
||
paModel,hasMemory,paCharacteristics,coefficients)
|
||
% 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 = ...
|
||
HelperNRNTNThroughput.getDefaultCoefficients;
|
||
paInputScaleFactor = -35;
|
||
else
|
||
hpa.CoefficientMatrix = coefficients;
|
||
end
|
||
hpaDelay = size(hpa.CoefficientMatrix,1)-1;
|
||
else
|
||
if paModel == "Custom"
|
||
if isempty(paCharacteristics)
|
||
hpa = comm.MemorylessNonlinearity(Method="Lookup table", ...
|
||
Table=HelperNRNTNThroughput.getDefaultLookup);
|
||
paInputScaleFactor = -35;
|
||
else
|
||
hpa = comm.MemorylessNonlinearity(Method="Lookup table", ...
|
||
Table=paCharacteristics);
|
||
end
|
||
elseif paModel == "2.1GHz GaAs"
|
||
hpa = @(in) HelperNRNTNThroughput.paMemorylessGaAs2Dot1GHz(in);
|
||
elseif paModel == "2.1GHz GaN"
|
||
hpa = @(in) HelperNRNTNThroughput.paMemorylessGaN2Dot1GHz(in);
|
||
elseif paModel == "28GHz CMOS"
|
||
hpa = @(in) HelperNRNTNThroughput.paMemorylessCMOS28GHz(in);
|
||
else % "28GHz GaN"
|
||
hpa = @(in) HelperNRNTNThroughput.paMemorylessGaN28GHz(in);
|
||
end
|
||
end
|
||
end
|
||
|
||
function params = computeAccessSlotParameters(simParameters,waveformInfo)
|
||
%computeAccessSlotParameters Finds the earliest access time
|
||
%between a moving UE and a satellite, and computes elevation
|
||
%angle, path loss, Doppler shift, and propagation delay for
|
||
%each slot.
|
||
|
||
% Initialize scenario sample time and stop time
|
||
ScenarioStopTime = simParameters.ScenarioStartTime + seconds(simParameters.ScenarioRunTime);
|
||
tSlot = HelperNRNTNThroughput.getSlotTimes( ...
|
||
simParameters.Carrier.SymbolsPerSlot,waveformInfo.SymbolLengths, ...
|
||
waveformInfo.SampleRate,simParameters.NFrames,simParameters.DataType);
|
||
ScenarioSampleTime = min(diff(tSlot));
|
||
% Create a satelliteScenario object
|
||
sc = satelliteScenario(simParameters.ScenarioStartTime,ScenarioStopTime,double(ScenarioSampleTime));
|
||
|
||
% Add a moving UE to the scenario
|
||
ueStartLLA = itrf2geographic(simParameters.UEPosition);
|
||
ueEndPos = simParameters.UEPosition+simParameters.UEVelocity*simParameters.ScenarioRunTime;
|
||
ueEndLLA = itrf2geographic(ueEndPos);
|
||
trajectory = geoTrajectory([ueStartLLA'; ueEndLLA'],[0 simParameters.ScenarioRunTime]);
|
||
gs = platform(sc,trajectory,"Name","UE");
|
||
|
||
% Add a satellite to the scenario using Keplerian elements
|
||
sat = satellite(sc,simParameters.SemiMajorAxis, ...
|
||
simParameters.Eccentricity, simParameters.Inclination, ...
|
||
simParameters.RAAN, simParameters.Argofperiapsis, ...
|
||
simParameters.TrueAnomaly, Name = "Satellite 1");
|
||
|
||
% Compute required access duration to transmit the entire data
|
||
totalTime = simParameters.NFrames*10e-3; % 10msec per frame
|
||
% Find the Access intervals
|
||
ac = access(sat, gs);
|
||
accessIntrvl = accessIntervals(ac);
|
||
if isempty(accessIntrvl)
|
||
error("No access intervals found: The satellite is not visible to the UE. " + ...
|
||
"Check satellite orbital parameters, UE position and velocity, " + ...
|
||
"and scenario start time and duration.");
|
||
end
|
||
% Find first interval with sufficient duration
|
||
disp('Access interval table (periods when satellite is visible to UE):');
|
||
disp(accessIntrvl)
|
||
|
||
if simParameters.LinkVariationModel=="Static"
|
||
accessStartTime = accessIntrvl.StartTime(1);
|
||
else
|
||
accessStartTime = [];
|
||
for idx=1:size(accessIntrvl,2)
|
||
if totalTime<=seconds(accessIntrvl.EndTime(idx)-accessIntrvl.StartTime(idx))
|
||
accessStartTime = accessIntrvl.StartTime(idx);
|
||
break;
|
||
end
|
||
end
|
||
|
||
if isempty(accessStartTime)
|
||
error("Insufficient visible time: The satellite's visibility duration is not long enough for continuous data transmission.")
|
||
end
|
||
end
|
||
disp("Transmission start time : " +string(accessStartTime) + " UTC")
|
||
|
||
% Calculate slot times
|
||
if simParameters.LinkVariationModel ~= "Time-varying"
|
||
slotTimes = accessStartTime;
|
||
else
|
||
tSlot = HelperNRNTNThroughput.getSlotTimes( ...
|
||
simParameters.Carrier.SymbolsPerSlot, ...
|
||
waveformInfo.SymbolLengths, ...
|
||
waveformInfo.SampleRate, ...
|
||
simParameters.NFrames, simParameters.DataType);
|
||
|
||
slotTimes = accessStartTime + seconds(tSlot);
|
||
end
|
||
|
||
% Preallocate
|
||
nSlots = numel(slotTimes);
|
||
el = zeros(nSlots, 1);
|
||
pathLoss = zeros(nSlots, 1);
|
||
preCompDoppler = zeros(nSlots, 1);
|
||
satDoppler = zeros(nSlots, 1);
|
||
delayInSeconds = zeros(nSlots, 1);
|
||
|
||
% Wavelength
|
||
c = physconst("lightspeed");
|
||
lambda = c / simParameters.CarrierFrequency;
|
||
uevel = [0;0;0]; % Stationary UE
|
||
|
||
for k = 1:nSlots
|
||
% Elevation angle
|
||
[~, el(k)] = aer(gs, sat, slotTimes(k));
|
||
% Doppler due to satellite and UE motion
|
||
preCompDoppler(k) = simParameters.EnableDoppler*dopplershift(gs, sat, slotTimes(k), ...
|
||
Frequency=simParameters.CarrierFrequency);
|
||
% Doppler shift only due to satellite motion
|
||
[satpos,satvel] = states(sat,slotTimes(k),CoordinateFrame='ecef');
|
||
uepos = states(gs,slotTimes(k),CoordinateFrame='ecef');
|
||
satDoppler(k) = simParameters.EnableDoppler*satcom.internal.dopplerShift(simParameters.CarrierFrequency,uepos,satpos,uevel,satvel);
|
||
% Delay
|
||
delay = latency(gs, sat, slotTimes(k));
|
||
delayInSeconds(k) = simParameters.EnableDelay*delay;
|
||
% Free space path loss
|
||
range = delay * c;
|
||
pathLoss(k) = simParameters.EnablePathLoss*fspl(double(range), lambda);
|
||
end
|
||
|
||
% Total number of slots
|
||
totalSlots = simParameters.NFrames*simParameters.Carrier.SlotsPerFrame;
|
||
|
||
% Expand to totalSlots if needed
|
||
if nSlots == 1
|
||
el = repmat(el, totalSlots, 1);
|
||
pathLoss = repmat(pathLoss, totalSlots, 1);
|
||
preCompDoppler = repmat(preCompDoppler, totalSlots, 1);
|
||
satDoppler = repmat(satDoppler, totalSlots, 1);
|
||
end
|
||
|
||
params.ElevationAngle = el;
|
||
params.PathLoss = pathLoss;
|
||
params.PreCompDopplerShift = preCompDoppler;
|
||
params.SatelliteDopplerShift = satDoppler;
|
||
params.DelayInSeconds = delayInSeconds;
|
||
params.SlotTimes = slotTimes;
|
||
params.AccessStartTime = accessStartTime;
|
||
end
|
||
|
||
function [txWaveform,info] = generatePUSCHWaveform(carrier,pusch,ulsch,wtx,dmrsPower,dt,numSlots)
|
||
%generatePUSCHWaveform generates a PUSCH (Physical Uplink
|
||
%Shared Channel) transmit waveform
|
||
|
||
arguments
|
||
carrier
|
||
pusch
|
||
ulsch = struct
|
||
wtx = eye(pusch.NumLayers)
|
||
dmrsPower = 3 % dB
|
||
dt = "double"
|
||
numSlots = 1
|
||
end
|
||
|
||
% Initialize variables
|
||
nTx = size(wtx,2); % Number of transmit antennas
|
||
tmpGrid = nrResourceGrid(carrier,nTx,OutputDataType=dt);
|
||
puschGrid = repmat(tmpGrid,[1 numSlots 1]);
|
||
refGrid = puschGrid;
|
||
nSlotSymb = carrier.SymbolsPerSlot;
|
||
initialNSlot = carrier.NSlot;
|
||
numCW = pusch.NumCodewords;
|
||
|
||
% Process loop for each slot
|
||
for slotIdx = 0:numSlots-1
|
||
[slotGrid,refSlotGrid] = deal(tmpGrid);
|
||
carrier.NSlot = initialNSlot + slotIdx;
|
||
|
||
% Generate PUSCH indices and PT-RS indices
|
||
[puschIndices,puschIndicesInfo,ptrsIndices] = nrPUSCHIndices(carrier,pusch);
|
||
if isempty(fieldnames(ulsch))
|
||
% Update the codeword with valid bits
|
||
cw = cell(1,numCW);
|
||
for i = 1:numCW
|
||
cw{i} = randi([0 1],puschIndicesInfo.G(i),1);
|
||
end
|
||
else
|
||
% Perform UL-SCH encoding with the inputs provided in the ulsch structure
|
||
cw = ulsch.Encoder(pusch.Modulation,pusch.NumLayers, ...
|
||
puschIndicesInfo.G,ulsch.RedundancyVersion, ...
|
||
ulsch.HARQProcessID,ulsch.CBGTI);
|
||
end
|
||
|
||
% Perform PUSCH modulation
|
||
if ~isempty(cw)
|
||
[puschSymbols, ptrsSymbols] = nrPUSCH(carrier,pusch,cw,OutputDataType=dt);
|
||
% Implementation-specific PUSCH MIMO precoding and mapping. This
|
||
% MIMO precoding step is in addition to any codebook based
|
||
% MIMO precoding done during PUSCH modulation above
|
||
[~,puschAntIndices] = nrExtractResources(puschIndices,slotGrid);
|
||
slotGrid(puschAntIndices) = puschSymbols * wtx;
|
||
end
|
||
|
||
% Implementation-specific PUSCH DM-RS MIMO precoding and mapping
|
||
% using the nrPDSCHPrecode function which supports PUSCH MIMO
|
||
% precoding as well. The first DM-RS creation includes codebook
|
||
% based MIMO precoding if applicable
|
||
dmrsSymbols = db2mag(dmrsPower)*nrPUSCHDMRS(carrier,pusch,OutputDataType=dt);
|
||
dmrsIndices = nrPUSCHDMRSIndices(carrier,pusch);
|
||
[dmrsAntSymbols,dmrsAntIndices] = nrPDSCHPrecode( ...
|
||
carrier,dmrsSymbols,dmrsIndices,wtx);
|
||
slotGrid(dmrsAntIndices) = dmrsAntSymbols;
|
||
refSlotGrid(dmrsAntIndices) = dmrsAntSymbols;
|
||
|
||
% Implementation-specific PUSCH PT-RS MIMO precoding and mapping
|
||
% using the nrPDSCHPrecode function which supports PUSCH MIMO
|
||
% precoding as well. If TransformPrecoding = 1, puschSymbols
|
||
% already include PT-RS symbols
|
||
if ~pusch.TransformPrecoding && ~isempty(ptrsSymbols)
|
||
[ptrsAntSymbols,ptrsAntIndices] = nrPDSCHPrecode(carrier,ptrsSymbols,ptrsIndices,wtx);
|
||
slotGrid(ptrsAntIndices) = ptrsAntSymbols;
|
||
refSlotGrid(ptrsAntIndices) = ptrsAntSymbols;
|
||
end
|
||
|
||
% Map to the grid containing multiple slots
|
||
symIdx = (nSlotSymb*slotIdx)+1:(nSlotSymb*(slotIdx+1));
|
||
puschGrid(:,symIdx,:) = slotGrid;
|
||
refGrid(:,symIdx,:) = refSlotGrid;
|
||
end
|
||
|
||
% Perform OFDM modulation
|
||
[txWaveform,info] = nrOFDMModulate(carrier,puschGrid);
|
||
|
||
% Append the resource grids to output structure
|
||
info.ResourceGrid = puschGrid;
|
||
info.ReferenceGrid = refGrid;
|
||
end
|
||
|
||
function [channel, maxChDelay, ueDoppler] = configureNTNChannel(simParameters,elevation,sampleRate)
|
||
%configureNTNChannel Configure NTN channel object and compute max channel delay.
|
||
|
||
mobileSpeed = norm(simParameters.UEVelocity); % m/s
|
||
c = physconst("LightSpeed");
|
||
ueDoppler = mobileSpeed*simParameters.CarrierFrequency/c;
|
||
|
||
if strcmpi(simParameters.NTNChannelType, 'Narrowband')
|
||
channel = p681LMSChannel;
|
||
channel.Environment = simParameters.Environment;
|
||
channel.AzimuthOrientation = simParameters.AzimuthOrientation;
|
||
channel.CarrierFrequency = simParameters.CarrierFrequency;
|
||
channel.ElevationAngle = elevation(1);
|
||
channel.MobileSpeed = mobileSpeed;
|
||
|
||
elseif strcmpi(simParameters.NTNChannelType, 'TDL')
|
||
channel = nrTDLChannel;
|
||
swapTransmitAndReceive(channel);
|
||
channel.TransmissionDirection = "Uplink";
|
||
channel.DelayProfile = simParameters.DelayProfile;
|
||
channel.DelaySpread = simParameters.DelaySpread;
|
||
channel.MaximumDopplerShift = ueDoppler;
|
||
channel.NumTransmitAntennas = simParameters.NumTransmitAntennas;
|
||
channel.NumReceiveAntennas = simParameters.NumReceiveAntennas;
|
||
|
||
elseif strcmpi(simParameters.NTNChannelType, 'CDL')
|
||
channel = nrCDLChannel;
|
||
swapTransmitAndReceive(channel);
|
||
channel.DelayProfile = simParameters.DelayProfile;
|
||
channel.SatelliteElevationAngle = elevation(1);
|
||
channel.DelaySpread = simParameters.DelaySpread;
|
||
channel.TransmitAntennaArray.Size = simParameters.TxArraySize;
|
||
channel.ReceiveAntennaArray.Size = simParameters.RxArraySize;
|
||
channel.CarrierFrequency = simParameters.CarrierFrequency;
|
||
|
||
else
|
||
error('Unknown NTNChannelType: %s', simParameters.NTNChannelType);
|
||
end
|
||
|
||
% Common channel parameters
|
||
channel.SampleRate = sampleRate;
|
||
channel.RandomStream = "mt19937ar with seed";
|
||
channel.Seed = 73;
|
||
|
||
% Maximum channel delay computation
|
||
chInfo = info(channel);
|
||
maxChDelay = ceil(max(chInfo.PathDelays * channel.SampleRate)) + ...
|
||
chInfo.ChannelFilterDelay;
|
||
|
||
% CDL-specific Doppler and orientation settings
|
||
if simParameters.NTNChannelType == "CDL"
|
||
channel.UTDirectionOfTravel = [chInfo.AnglesAoA(1) chInfo.AnglesAoD(1); ... % [Rx_Sat_Azi Tx_UE_Azi; Rx_Sat_Zen Tx_UE_Zen]
|
||
chInfo.AnglesZoA(1) chInfo.AnglesZoD(1)];
|
||
channel.ReceiveArrayOrientation = [chInfo.AnglesAoA(1); chInfo.AnglesZoA(1) - 90; 0]; % Satellite RX array pointing toward UE
|
||
end
|
||
end
|
||
|
||
function initialOffset = estimateInitialTimingOffset(rxCarrier,rxData,dmrsParams,syncModel,inifValsVec)
|
||
%estimateInitialTimingOffset Estimate initial timing offset based on sync model
|
||
if syncModel == "auto corr"
|
||
initialOffset = nrTimingEstimate(rxCarrier,rxData,dmrsParams.DMRSIndices,dmrsParams.DMRSSymbols);
|
||
elseif syncModel == "diff corr"
|
||
initialOffset = HelperNRNTNThroughput.diffcorr( ...
|
||
rxCarrier,rxData,dmrsParams.DMRSIndices,dmrsParams.DMRSSymbols);
|
||
else
|
||
initialOffset = HelperNRNTNThroughput.jointTimeFreq( ...
|
||
rxCarrier,rxData,dmrsParams.DMRSIndices,dmrsParams.DMRSSymbols,inifValsVec);
|
||
end
|
||
end
|
||
|
||
function [rxWaveform,offset,shiftOut] = compensateRxDoppler(rxWaveform,rxDopplerParams,dmrsParams, waveinfo)
|
||
%compensateRxDoppler Compensate for Doppler shift in received waveform
|
||
|
||
if rxDopplerParams.CompensationMethod == "joint time-freq"
|
||
% Perform joint time-frequency synchronization
|
||
[offset,fOEst] = HelperNRNTNThroughput.jointTimeFreq( ...
|
||
rxDopplerParams.RxCarrier,rxWaveform,dmrsParams.DMRSIndices,dmrsParams.DMRSSymbols,rxDopplerParams.FreqSearchRange);
|
||
% Compensate Doppler shift
|
||
rxWaveform = HelperNRNTNThroughput.compensateDopplerShift( ...
|
||
rxWaveform,waveinfo.SampleRate, ...
|
||
fOEst,true);
|
||
% Estimate and compensate the residual (fractional) Doppler shift
|
||
fractionalDopplerShift = ...
|
||
HelperNRNTNThroughput.estimateFractionalDopplerShift( ...
|
||
rxWaveform,rxDopplerParams.RxCarrier.SubcarrierSpacing,waveinfo.Nfft, ...
|
||
waveinfo.CyclicPrefixLengths(2),0,true);
|
||
rxWaveform = HelperNRNTNThroughput.compensateDopplerShift( ...
|
||
rxWaveform,waveinfo.SampleRate, ...
|
||
fractionalDopplerShift,true);
|
||
shiftOut = [];
|
||
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] = ...
|
||
HelperNRNTNThroughput.estimateFractionalDopplerShift( ...
|
||
rxWaveform,rxDopplerParams.RxCarrier.SubcarrierSpacing,waveinfo.Nfft, ...
|
||
waveinfo.CyclicPrefixLengths(2),rxDopplerParams.Threshold, ...
|
||
true);
|
||
rxWaveform = HelperNRNTNThroughput.compensateDopplerShift( ...
|
||
rxWaveform,waveinfo.SampleRate, ...
|
||
fractionalDopplerShift,true);
|
||
|
||
% Perform integer Doppler frequency shift estimation and
|
||
% compensation. Use the demodulation reference signals to
|
||
% compute the integer Doppler shift.
|
||
[integerDopplerShift,shiftOut] = ...
|
||
HelperNRNTNThroughput.estimateIntegerDopplerShift( ...
|
||
rxDopplerParams.RxCarrier,rxWaveform,dmrsParams.DMRSIndices, ...
|
||
dmrsParams.DMRSSymbols,rxDopplerParams.SampleOffset, ...
|
||
rxDopplerParams.UsePrevShift,rxDopplerParams.UseDiffCorrFlag,rxDopplerParams.Shift,rxDopplerParams.TotalDelay, ...
|
||
detFlag);
|
||
rxWaveform = HelperNRNTNThroughput.compensateDopplerShift( ...
|
||
rxWaveform,waveinfo.SampleRate, ...
|
||
integerDopplerShift,true);
|
||
offset = [];
|
||
end
|
||
end
|
||
|
||
function [decbits,blkerr,cbgerr,info] = recoverPUSCHBits(rxWaveform,rxCarrier,dmrsParams,pusch,ulschDecParams)
|
||
%recoverPUSCHBits Recovers bits from the received waveform
|
||
|
||
% 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 PUSCH DM-RS for each
|
||
% layer. This channel estimate includes the effect of transmitter
|
||
% precoding.
|
||
[estChannelGrid,noiseEst] = nrChannelEstimate(rxCarrier,rxGrid,...
|
||
dmrsParams.DMRSIndices,dmrsParams.DMRSSymbols,'CDMLengths',pusch.DMRS.CDMLengths);
|
||
|
||
[refPUSCHIndices,refPUSCHIndicesInfo] = nrPUSCHIndices(rxCarrier,pusch);
|
||
|
||
% Get PUSCH REs from the received grid and estimated channel grid
|
||
[puschRx,puschHest] = nrExtractResources(...
|
||
refPUSCHIndices,rxGrid,estChannelGrid);
|
||
|
||
% Perform equalization
|
||
[puschEq,csi] = nrEqualizeMMSE(puschRx,puschHest,noiseEst);
|
||
|
||
% Common phase error (CPE) compensation
|
||
if pusch.EnablePTRS
|
||
puschEq = hCompensateCPE(rxCarrier,pusch,puschEq,rxGrid,estChannelGrid,noiseEst);
|
||
end
|
||
|
||
if ~any(isnan(puschEq))
|
||
% Decode PUSCH symbols
|
||
[ulschLLRs,rxSymbols] = nrPUSCHDecode(rxCarrier,pusch,puschEq,noiseEst);
|
||
|
||
% Apply channel state information to LLRs supporting transform precoding & PT-RS.
|
||
ulschLLRs = HelperNRNTNThroughput.applyCSIToLLRs(rxCarrier, pusch, numel(refPUSCHIndicesInfo.PRBSet), rxSymbols, ulschLLRs, csi);
|
||
|
||
% Decode the UL-SCH transport channel
|
||
[decbits,blkerr,cbgerr] = ulschDecParams.Decoder(ulschLLRs,pusch.Modulation,...
|
||
pusch.NumLayers,ulschDecParams.RedundancyVersion,...
|
||
ulschDecParams.HARQProcessID,ulschDecParams.CBGTI);
|
||
else
|
||
decbits = [];
|
||
blkerr = true;
|
||
cbgerr = true;
|
||
end
|
||
|
||
info = refPUSCHIndicesInfo;
|
||
end
|
||
|
||
function ulschLLRs = applyCSIToLLRs(rxCarrier, pusch, MRB, rxSymbols, ulschLLRs, csi)
|
||
%applyCSIToLLRs Apply channel state information to LLRs,
|
||
%supporting transform precoding & PT-RS.
|
||
|
||
% Ensure LLRs and symbols are cell arrays (one per codeword)
|
||
if ~iscell(ulschLLRs)
|
||
ulschLLRs = num2cell(ulschLLRs, 1);
|
||
rxSymbols = num2cell(rxSymbols, 1);
|
||
end
|
||
|
||
% Apply channel state information (CSI) produced by the equalizer,
|
||
% including the effect of transform precoding if enabled
|
||
if pusch.TransformPrecoding
|
||
MSC = MRB * 12; % Number of subcarriers
|
||
|
||
% Undo transform precoding for CSI, normalize
|
||
csi = nrTransformDeprecode(csi, MRB) / sqrt(MSC);
|
||
% Repeat each de-precoded CSI value across its corresponding subcarriers
|
||
csi = repmat(csi((1:MSC:end).'), 1, MSC).';
|
||
if pusch.EnablePTRS
|
||
% Remove PT-RS CSI
|
||
ptrsLayerIndices = nrPUSCHPTRSIndices(rxCarrier, pusch);
|
||
csi(ptrsLayerIndices) = [];
|
||
end
|
||
|
||
% Reshape CSI to match the shape of rxSymbols
|
||
csi = reshape(csi, size(rxSymbols{:}));
|
||
end
|
||
|
||
% Layer de-map CSI
|
||
csi = nrLayerDemap(csi);
|
||
|
||
% Apply CSI to LLRs for each codeword
|
||
for cwIdx = 1:pusch.NumCodewords
|
||
Qm = length(ulschLLRs{cwIdx}) / length(rxSymbols{cwIdx}); % Bits per symbol
|
||
csi{cwIdx} = repmat(csi{cwIdx}.', Qm, 1); % Expand by each bit per symbol
|
||
ulschLLRs{cwIdx} = ulschLLRs{cwIdx} .* csi{cwIdx}(:); % Scale by CSI
|
||
end
|
||
end
|
||
|
||
end
|
||
|
||
end
|
||
|
||
function lla = itrf2geographic(ecef)
|
||
%itrf2geographic converts Earth-Centered Earth-Fixed (ECEF)
|
||
%coordinates to geodetic coordinates: latitude (degrees),
|
||
%longitude (degrees), and altitude (meters).
|
||
|
||
geographicCoordinates = matlabshared.orbit.internal.Transforms.itrf2geographic(ecef);
|
||
% Convert to degrees and get the LLA coordinates
|
||
lla = [rad2deg(geographicCoordinates(1,1:end)); ...
|
||
rad2deg(geographicCoordinates(2,1:end)); ...
|
||
geographicCoordinates(3,1:end)];
|
||
end
|