diff --git a/lab_5/HelperNRNTNThroughput.m b/lab_5/HelperNRNTNThroughput.m new file mode 100755 index 0000000..1859ed5 --- /dev/null +++ b/lab_5/HelperNRNTNThroughput.m @@ -0,0 +1,1221 @@ +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 + 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 diff --git a/lab_5/HelperNTNPositionEstimate.m b/lab_5/HelperNTNPositionEstimate.m new file mode 100755 index 0000000..71911b1 --- /dev/null +++ b/lab_5/HelperNTNPositionEstimate.m @@ -0,0 +1,151 @@ +classdef HelperNTNPositionEstimate + %HelperNTNPositionEstimate Class defining the supporting functions used + %in the receiver position estimation in NTN example + % + % Note: This is an undocumented class and its API and/or + % functionality may change in subsequent releases. + + % Copyright 2024-2025 The MathWorks, Inc. + + methods (Static) + function out = instantposition(s, v, r, rDot) + % OUT = instantposition(S,V,R,RDOT) returns the possible + % position OUT using satellite position S, satellite velocity + % V, range R, and range-rate RDOT. S and V are in Earth + % Centered Earth Fixed (ECEF) coordinate system. R is in m and + % RDOT is in m/s. OUT is a 3-by-2 matrix with each column + % corresponding to a possible position in LLA (Latitude, + % Longitude, Altitude). + + tolerance = 1e-4; + max_iter_refine = 1000; + % Earth Parameters + r_E = 6378137; % Equatorial radius (meters) + r_p = 6356752; % Polar radius (meters) + + % Initialize two r_c values with random values + r_c_1 = (r_E - r_p) * rand(); + r_c_2 = r_c_1; % Starting both from the same value (can be different if needed) + + % Refinement loop for two solutions + for iter = 1:max_iter_refine + prev_r_c_1 = r_c_1; + prev_r_c_2 = r_c_2; + + % Compute geometric parameters + r_s2 = dot(s,s); % Squared norm of position vector + + % Linear system for refinement + A = [s(1) s(2); v(1) v(2)]; + b_1 = 0.5*(r_s2 + r_c_1^2 - r^2); + b_2 = 0.5*(r_s2 + r_c_2^2 - r^2); + c = dot(s,v) - r*rDot; + + % Solve the parameters for r_c_1 and r_c_2 + mat_1 = [b_1, s(3); c, v(3)]; + mat_2 = [b_2, s(3); c, v(3)]; + + params_1 = A \ mat_1; + params_2 = A \ mat_2; + + % Extract parameters for the first solution + alpha_1 = params_1(1); + beta_1 = params_1(3); + gamma_1 = params_1(2); + delta_1 = params_1(4); + + % Extract parameters for the second solution + alpha_2 = params_2(1); + beta_2 = params_2(3); + gamma_2 = params_2(2); + delta_2 = params_2(4); + + % Calculate discriminant for both solutions + numerator_1 = alpha_1 * beta_1 + gamma_1 * delta_1; + discriminant_1 = numerator_1^2 - (1 + delta_1^2 + beta_1^2) * (alpha_1^2 + gamma_1^2 - r_c_1^2); + + numerator_2 = alpha_2 * beta_2 + gamma_2 * delta_2; + discriminant_2 = numerator_2^2 - (1 + delta_2^2 + beta_2^2) * (alpha_2^2 + gamma_2^2 - r_c_2^2); + + % Handle negative discriminant by regenerating r_c + if discriminant_1 < 0 + r_c_1 = r_p + (r_E - r_p) * rand(); + continue; + end + if discriminant_2 < 0 + r_c_2 = r_p + (r_E - r_p) * rand(); + continue; + end + + % Compute z-coordinates for both solutions + z_hat_1 = (numerator_1 + sqrt(discriminant_1)) / (1 + delta_1^2 + beta_1^2); + z_hat_2 = (numerator_2 - sqrt(discriminant_2)) / (1 + delta_2^2 + beta_2^2); + + % Compute x and y coordinates for both solutions + x_hat_1 = alpha_1 - beta_1 * z_hat_1; + y_hat_1 = gamma_1 - delta_1 * z_hat_1; + + x_hat_2 = alpha_2 - beta_2 * z_hat_2; + y_hat_2 = gamma_2 - delta_2 * z_hat_2; + + % Convert to LLA + lla_1 = ecef2llaLocal([x_hat_1; y_hat_1; z_hat_1]); + lla_2 = ecef2llaLocal([x_hat_2; y_hat_2; z_hat_2]); + + % Recalculate r_c for refined solutions + dr_1 = sind(lla_1(1))^2 + cosd(lla_1(1))^2 * (r_p / r_E)^2; + r_c_1 = sqrt(r_p^2 / dr_1); + + dr_2 = sind(lla_2(1))^2 + cosd(lla_2(1))^2 * (r_p / r_E)^2; + r_c_2 = sqrt(r_p^2 / dr_2); + + % Check convergence for both solutions + if abs(r_c_1 - prev_r_c_1) < tolerance && abs(r_c_2 - prev_r_c_2) < tolerance + break; + end + end + out = []; + if discriminant_1 > 0 + out = lla_1(:); + end + if discriminant_2 > 0 + out = [out(:) lla_2(:)]; + end + end + + function error = positionerror(originalPos,estimatedPos) + % ERR = positionerror(ORIGINALPOS,ESTIMATEDPOS) computes the + % position errors ERR for given original position ORIGINALPOS + % and estimated positions ESTIMATEDPOS. Original position and + % estimated positions are in LLA. ESTIMATEDPOS is 3-by-N matrix + % with N being the number of positions. + + % Convert the positions to ECEF + posECEF = lla2ecefLocal([originalPos estimatedPos]); + error = vecnorm(posECEF(:,1)-posECEF(:,2:end)); + end + end +end + +function lla = ecef2llaLocal(ecef) +% Output Latitude (degrees), Longitude (degrees), and Altitude (meters) + +tmp = matlabshared.orbit.internal.Transforms.itrf2geographic(ecef); +% Convert to degrees and get the LLA coordinates +lla = [rad2deg(tmp(1,1:end)); ... + rad2deg(tmp(2,1:end)); ... + tmp(3,1:end)]; + +end + +function ecef = lla2ecefLocal(lla) +% Output ECEF: X (meters), Y (meters), and Z (meters) + +% Convert to radians +lla = [deg2rad(lla(1,1:end)); ... + deg2rad(lla(2,1:end)); ... + lla(3,1:end)]; +% Get ECEF coordinates +ecef = matlabshared.orbit.internal.Transforms.geographic2itrf(lla); + +end diff --git a/lab_5/create_baseline_configuration.m b/lab_5/create_baseline_configuration.m new file mode 100644 index 0000000..1bcd69a --- /dev/null +++ b/lab_5/create_baseline_configuration.m @@ -0,0 +1,73 @@ +function [cfgEHT, simParameters, channel, chInfo, maxChDelay, satelliteDopplerShift, chanBW] = create_baseline_configuration(profile, ... + bandwidth_mhz, ... + include_free_space_pathloss, ... + delay_model, ... + carrier_frequency, ... + mobile_speed, ... + mobile_altitude, ... + satellite_altitude, ... + random_stream, ... + seed, ... + tx_antenna_count, ... + rx_antenna_count, ... + apep) + + switch bandwidth_mhz + case 20 + chanBW = 'CBW20'; + case 40 + chanBW = 'CBW40'; + case 80 + chanBW = 'CBW80'; + case 160 + chanBW = 'CBW160'; + case 320 + chanBW = 'CBW320'; + otherwise + error("Not supported bandwidth!"); + end + + cfgEHT = wlanEHTMUConfig(chanBW); + cfgEHT.User{1}.APEPLength = apep; + cfgEHT.NumTransmitAntennas = tx_antenna_count; + cfgEHT.User{1}.NumSpaceTimeStreams = 1; + cfgEHT.User{1}.MCS = 0; % default to sweep later + fs = wlanSampleRate(chanBW); + + simParameters.IncludeFreeSpacePathLoss = include_free_space_pathloss; + simParameters.DelayModel = delay_model; + simParameters.CarrierFrequency = carrier_frequency; + simParameters.ElevationAngle = 0; % default to sweep later + simParameters.MobileSpeed = mobile_speed; + simParameters.MobileAltitude = mobile_altitude; + simParameters.SatelliteAltitude = satellite_altitude; + simParameters.SampleRate = fs; + simParameters.RandomStream = random_stream; + simParameters.Seed = seed; + simParameters.NumTransmitAntennas = tx_antenna_count; + simParameters.NumReceiveAntennas = rx_antenna_count; + + simParameters.NTNChannelType = "TDL"; + simParameters.DelayProfile = "NTN-TDL-" + profile; + simParameters.DelaySpread = 30e-9; + + c = physconst("lightspeed"); + satelliteDopplerShift = dopplerShiftCircularOrbit( ... + simParameters.ElevationAngle,simParameters.SatelliteAltitude, ... + simParameters.MobileAltitude,simParameters.CarrierFrequency); + + 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; + + channel.SampleRate = simParameters.SampleRate; + channel.RandomStream = simParameters.RandomStream; + channel.Seed = simParameters.Seed; + + chInfo = info(channel); + maxChDelay = ceil(max(chInfo.PathDelays*channel.SampleRate)) + chInfo.ChannelFilterDelay; +end \ No newline at end of file diff --git a/lab_5/create_satellite_sim.m b/lab_5/create_satellite_sim.m new file mode 100644 index 0000000..82127c5 --- /dev/null +++ b/lab_5/create_satellite_sim.m @@ -0,0 +1,12 @@ +function satSim = create_satellite_sim(satConfig, stationConfig) + satSim.sc = satelliteScenario; + satSim.sc.StartTime = satConfig.startTime; + satSim.sc.StopTime = satConfig.stopTime; + satSim.sc.SampleTime = satConfig.sampleTime; + + earthRadius = physconst("EarthRadius"); + semiMajorAxis = earthRadius + satConfig.satAltitude; + + satSim.sat = satellite(satSim.sc, semiMajorAxis, satConfig.eccentricity, satConfig.inclination, satConfig.raan, satConfig.argOfPeriapsis, satConfig.trueAnomaly, Name=satConfig.name, OrbitPropagator=satConfig.orbitPropagator); + satSim.ue = groundStation(satSim.sc, stationConfig.Lat, stationConfig.Lon, Altitude=stationConfig.Altitude, Name=stationConfig.Name); +end \ No newline at end of file diff --git a/lab_5/distance_over_time.svg b/lab_5/distance_over_time.svg new file mode 100644 index 0000000..dca0248 --- /dev/null +++ b/lab_5/distance_over_time.svg @@ -0,0 +1,207 @@ + + +Qt SVG Document +MATLAB, The MathWorks, Inc. Version 26.1.0.3251617 R2026a Update 2 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/lab_5/elevation_angle_over_time.svg b/lab_5/elevation_angle_over_time.svg new file mode 100644 index 0000000..b5d014f --- /dev/null +++ b/lab_5/elevation_angle_over_time.svg @@ -0,0 +1,226 @@ + + +Qt SVG Document +MATLAB, The MathWorks, Inc. Version 26.1.0.3251617 R2026a Update 2 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +08:00 + + + +08:30 + + + +09:00 + + + +09:30 + + + +10:00 + + + +Time [s] + + + +Jun 17, 2026 + + + +Elevation angle to satellite [degree] + + + + + + + + + + + + + + + + + + + + + + + + + + + + +-100 + + + +-80 + + + +-60 + + + +-40 + + + +-20 + + + +0 + + + +20 + + + +40 + + + +60 + + + +80 + + + +Elevation angle over time plot of satellite pass + + + + + + + diff --git a/lab_5/lab_5.m b/lab_5/lab_5.m new file mode 100644 index 0000000..591cd3f --- /dev/null +++ b/lab_5/lab_5.m @@ -0,0 +1,31 @@ +satConfig.startTime = datetime(2026,6,17,08,00,00); +satConfig.stopTime = satConfig.startTime + hours(2); +satConfig.sampleTime = 30; +satConfig.satAltitude = 600e3; +satConfig.eccentricity = 0; +satConfig.inclination = 90; +satConfig.raan = 42; +satConfig.argOfPeriapsis = 32; +satConfig.trueAnomaly = 0; +satConfig.name = "Timos Surveillance LEO Satellite"; +satConfig.orbitPropagator = "two-body-keplerian"; + +groundStation.Lat = 52.515744683039436; +groundStation.Lon = 13.325825073341306; +groundStation.Name = "HFT building at Einsteinufer 25"; +% https://www.randymajors.org/elevation-on-google-maps?x=13.3258251&y=52.5157447&cx=13.3258251&cy=52.5157447&zoom=15&labels=show +groundStation.Altitude = 37; + +calc_task1 = true; +calc_task2 = true; + +if calc_task1 + task1; +end + +if calc_task2 + if ~calc_task1 + task1; + end + task2; +end \ No newline at end of file diff --git a/lab_5/leo_satellite_pass.md b/lab_5/leo_satellite_pass.md new file mode 100644 index 0000000..3528c5e --- /dev/null +++ b/lab_5/leo_satellite_pass.md @@ -0,0 +1,29 @@ +

+Lab 5 - Wireless Networking Technologies +

+

+Analyzing a LEO Satellite Pass +

+

+by Timo Niemann +

+ +### Task1 + +#### 1.1 + +See code `lab_5.m` and `create_satellite_sim.m`, where the first file defines all parameters for the simulation and the create_satellite_sim function creates and returns the satellite and ground station. + +#### 1.2 + +To calculate the distance to the satellite, the MATLAB access analysis example was used as a basis and modified. The distance is calculated with the `aer` function, which returns azimuth, elevation, and slant range between the ground station and the satellite. To avoid invalid channel simulations, the signal-based estimation step is skipped when the Doppler shift `fshift` is non-finite. + +![Distance over time sweep](distance_over_time.svg) + +The first satellite flyby occurs at around 08:05:30. During the pass, the distance first decreases as the satellite approaches the ground station and then increases again after the closest approach. Around 08:53:30 to 08:54:00, the distance reaches a maximum because the satellite is on the far side of its orbit relative to the ground station. The second pass occurs at around 09:42. Its minimum distance is larger than during the first pass because Earth's rotation changes the relative position of the ground station below the satellite orbit. + +#### 1.3 + +![Elevation angle over time sweep](elevation_angle_over_time.svg) + +As in the distance plot, the elevation angle reaches its highest value during the first pass. The second pass has a lower maximum elevation because Earth's rotation shifts the ground station relative to the satellite ground track. diff --git a/lab_5/parfor_waitbar.m b/lab_5/parfor_waitbar.m new file mode 100644 index 0000000..5fa352c --- /dev/null +++ b/lab_5/parfor_waitbar.m @@ -0,0 +1,29 @@ +function parfor_waitbar(f, totalJobs, detail, reset) + persistent done + + if reset + done = 0; + if ishandle(f) + if isempty(detail) + waitbar(0, f, "Starting"); + else + waitbar(0, f, detail); + end + end + return + end + + if isempty(done) + done = 0; + end + + done = done + 1; + + if ishandle(f) + progressText = sprintf("Done %d/%d", done, totalJobs); + if ~isempty(detail) + progressText = sprintf("%s | %s", progressText, char(detail)); + end + waitbar(min(done / totalJobs, 1), f, progressText); + end +end diff --git a/lab_5/perform_access_analysis.m b/lab_5/perform_access_analysis.m new file mode 100644 index 0000000..825beeb --- /dev/null +++ b/lab_5/perform_access_analysis.m @@ -0,0 +1,75 @@ +function res = perform_access_analysis(satSim, queryTime) + sat = satSim.sat; + ue = satSim.ue; + + ac = access(sat,ue); + [status,timeout] = accessStatus(ac); + statusIdx = find(status,1); + if isempty(statusIdx) + disp("Position estimation cannot be performed " + ... + "as there is no access between satellite and UE.") + return + end + delete(ac) + + [~,el,r] = aer(ue,sat,queryTime); + + res.distanceToSatellite = r; + res.elevationAngleToSatellite = el; + + carrierFrequency = 2e9; + [fshift,~,dopplerInfo] = dopplershift( ... + sat,ue,queryTime,Frequency=carrierFrequency); + res.rangeRate = -1*dopplerInfo.RelativeVelocity; + + [satPos,satVel] = states(sat,queryTime,CoordinateFrame="ecef"); + + res.satPos = satPos; + res.satVel = satVel; + + if ~isfinite(fshift) + return; + end + + rng(0) + lightspeed = physconst("lightspeed"); + + frameDurInSeconds = 10e-3; + delayInSeconds = r/lightspeed; + nFrames = ceil(delayInSeconds/frameDurInSeconds) + 1; + carrier = nrCarrierConfig; + pdsch = nrPDSCHConfig; + pdsch.DMRS.DMRSAdditionalPosition = 2; + [txWaveform,waveformInfo] = ... + HelperNRNTNThroughput.generatePDSCHWaveform( ... + carrier,pdsch,struct,eye(pdsch.NumLayers),"double",carrier.SlotsPerFrame*nFrames); + refGridAllSlots = waveformInfo.ResourceGrid; + + delayInSamples = delayInSeconds.*waveformInfo.SampleRate; + integDelaySamples = floor(delayInSamples); + fracDelaySamples = (delayInSamples - integDelaySamples); + staticDelay = dsp.Delay(Length=integDelaySamples); + variableFractionalDelay = dsp.VariableFractionalDelay( ... + InterpolationMethod="Farrow", ... + FarrowSmallDelayAction="Use off-centered kernel", ... + MaximumDelay=1); + + delayedTx = staticDelay(txWaveform); + delayedTx = variableFractionalDelay(delayedTx,fracDelaySamples); + + channel = p681LMSChannel; + channel.SampleRate = waveformInfo.SampleRate; + channel.CarrierFrequency = carrierFrequency; + channel.SatelliteDopplerShift = fshift; + channel.MobileSpeed = 0; + channel.ElevationAngle = el; + channel.RandomStream = "mt19937ar with seed"; + rxWaveform = channel(delayedTx); + + freqSearch = 0.999*fshift:sign(fshift):1.001*fshift; + [delayEstSamples,fshiftEst] = HelperNRNTNThroughput.jointTimeFreq( ... + carrier,rxWaveform,refGridAllSlots,freqSearch); + + res.rangeEst = delayEstSamples*lightspeed/waveformInfo.SampleRate; + res.rangeRateEst = -1*fshiftEst*lightspeed/carrierFrequency; +end \ No newline at end of file diff --git a/lab_5/receive_power_over_time.svg b/lab_5/receive_power_over_time.svg new file mode 100644 index 0000000..e426239 --- /dev/null +++ b/lab_5/receive_power_over_time.svg @@ -0,0 +1,262 @@ + + +Qt SVG Document +MATLAB, The MathWorks, Inc. Version 26.1.0.3251617 R2026a Update 2 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +08:00 + + + +08:30 + + + +09:00 + + + +09:30 + + + +10:00 + + + +Time + + + +Jun 17, 2026 + + + +Power [dBm] + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +-170 + + + +-160 + + + +-150 + + + +-140 + + + +-130 + + + +-120 + + + +-110 + + + +-100 + + + +-90 + + + +-80 + + + +-70 + + + +Receive power over time plot of satellite pass + + + + + + + + + + + +receive sensitivity -85 dBm + + + + + + + +reliable reception -75 dBm + + + diff --git a/lab_5/required_transmit_power.svg b/lab_5/required_transmit_power.svg new file mode 100644 index 0000000..64554e9 --- /dev/null +++ b/lab_5/required_transmit_power.svg @@ -0,0 +1,196 @@ + + +Qt SVG Document +MATLAB, The MathWorks, Inc. Version 26.1.0.3251617 R2026a Update 2 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +08:00 + + + +08:30 + + + +09:00 + + + +09:30 + + + +10:00 + + + +Time + + + +Jun 17, 2026 + + + +Power [dBm] + + + + + + + + + + + + + + + + + + + + + + +80 + + + +85 + + + +90 + + + +95 + + + +100 + + + +105 + + + +110 + + + +Required transmit power for -75 dBm receive power + + + + + + + diff --git a/lab_5/simulateTransmission.m b/lab_5/simulateTransmission.m new file mode 100644 index 0000000..39734a4 --- /dev/null +++ b/lab_5/simulateTransmission.m @@ -0,0 +1,106 @@ +function [numPacketErrors, numPacketsSimulated, cfrMagnitudes, rxPower, pathLoss]= simulateTransmission(accessAnalysisResults, cfgEHT, channel, maxChDelay, snr, maxNumErrors, simParameters) + numPacketErrors = 0; + numPacketsSimulated = 0; + cfrMagnitudes = []; + rxPower = zeros(size(accessAnalysisResults.elevationAngles)); + pathLoss = zeros(size(accessAnalysisResults.elevationAngles)); + + ofdmInfo = wlanEHTOFDMInfo('EHT-Data', cfgEHT); + ind = wlanFieldIndices(cfgEHT); + + stream = RandStream('combRecursive', Seed=simParameters.Seed); + stream.Substream = 1; + RandStream.setGlobalStream(stream); + + powerOnlySimulation = isempty(snr); + if ~powerOnlySimulation + snrVal = convertSNR(snr, "snrsc", "snr", FFTLength=ofdmInfo.FFTLength, NumActiveSubcarriers=ofdmInfo.NumTones); + end + + if isfield(simParameters, 'TransmitPower') + txPower = simParameters.TransmitPower; + else + txPower = 20; + end + + lambda = physconst("lightspeed") / simParameters.CarrierFrequency; + includePathLoss = isfield(simParameters, 'IncludeFreeSpacePathLoss') && simParameters.IncludeFreeSpacePathLoss; + + while (powerOnlySimulation || numPacketErrors < maxNumErrors) && numPacketsSimulated < numel(accessAnalysisResults.elevationAngles) + numPacketsSimulated = numPacketsSimulated + 1; + + if includePathLoss + packetDistance = accessAnalysisResults.distances(numPacketsSimulated); + pathLoss(numPacketsSimulated) = 20 * log10(4 * pi * packetDistance / lambda); + else + pathLoss(numPacketsSimulated) = 0; + end + rxPower(numPacketsSimulated) = txPower - pathLoss(numPacketsSimulated); + + txPSDU = randi([0 1], psduLength(cfgEHT) * 8, 1); + tx = wlanWaveformGenerator(txPSDU, cfgEHT); + + stf = tx(ind.LSTF(1):ind.LSTF(2), :); + stfPower = mean(abs(stf(:)).^2); + + tx = tx * sqrt(10^((txPower - 30) / 10)) / sqrt(stfPower); + + txPad = [tx; zeros(maxChDelay, cfgEHT.NumTransmitAntennas)]; + + if powerOnlySimulation + continue; + end + + packetElevationAngle = accessAnalysisResults.elevationAngles(numPacketsSimulated); + channel.SatelliteDopplerShift = dopplerShiftCircularOrbit(packetElevationAngle, simParameters.SatelliteAltitude, simParameters.MobileAltitude, simParameters.CarrierFrequency); + + reset(channel); + rx = awgn(channel(txPad) * db2mag(-pathLoss(numPacketsSimulated)), snrVal); + + coarsePacketOffset = wlanPacketDetect(rx, cfgEHT.ChannelBandwidth); + if isempty(coarsePacketOffset) + numPacketErrors = numPacketErrors + 1; + continue; + end + + lstf = rx(coarsePacketOffset + (ind.LSTF(1) : ind.LSTF(2)), :); + coarseFrequencyOffset = wlanCoarseCFOEstimate(lstf, cfgEHT.ChannelBandwidth); + rx = frequencyOffset(rx, simParameters.SampleRate, -coarseFrequencyOffset); + + nonhtfields = rx(coarsePacketOffset + (ind.LSTF(1) : ind.LSIG(2)), :); + finePacketOffset = wlanSymbolTimingEstimate(nonhtfields, cfgEHT.ChannelBandwidth); + + packetOffset = coarsePacketOffset + finePacketOffset; + + if packetOffset > maxChDelay + numPacketErrors = numPacketErrors + 1; + continue; + end + + rxLLTF = rx(packetOffset + (ind.LLTF(1) : ind.LLTF(2)), :); + fineFrequencyOffset = wlanFineCFOEstimate(rxLLTF, cfgEHT.ChannelBandwidth); + rx = frequencyOffset(rx, simParameters.SampleRate, -fineFrequencyOffset); + + rxHELTF = rx(packetOffset + (ind.EHTLTF(1) : ind.EHTLTF(2)), :); + heltfDemod = wlanEHTDemodulate(rxHELTF, "EHT-LTF", cfgEHT); + [chanEstimate, pilotEstimate] = wlanEHTLTFChannelEstimate(heltfDemod, cfgEHT); + + cfrMagnitudes(:, end + 1) = abs(chanEstimate(:, 1, 1)); + + rxData = rx(packetOffset + (ind.EHTData(1) : ind.EHTData(2)), :); + demodSym = wlanEHTDemodulate(rxData, "EHT-Data", cfgEHT); + + demodSym = wlanEHTTrackPilotError(demodSym, chanEstimate, cfgEHT, "EHT-Data"); + + nVarEst = wlanEHTDataNoiseEstimate(demodSym(ofdmInfo.PilotIndices, :, :),pilotEstimate, cfgEHT); + + demodDataSym = demodSym(ofdmInfo.DataIndices, :, :); + chanEstimateData = chanEstimate(ofdmInfo.DataIndices, :, :); + + [eqSym, csi] = wlanEHTEqualize(demodDataSym, chanEstimateData, nVarEst, cfgEHT, 'EHT-Data'); + + rxPSDU = wlanEHTDataBitRecover(eqSym, nVarEst, csi, cfgEHT); + + numPacketErrors = numPacketErrors + any(biterr(txPSDU, rxPSDU)); + end +end diff --git a/lab_5/task1.m b/lab_5/task1.m new file mode 100644 index 0000000..b19f1d9 --- /dev/null +++ b/lab_5/task1.m @@ -0,0 +1,45 @@ +satSim = create_satellite_sim(satConfig, groundStation); +% viewer = satelliteScenarioViewer(satSim.sc); + +sampleTimes = satConfig.startTime:seconds(satConfig.sampleTime):satConfig.stopTime; +distancesToSat = zeros(size(sampleTimes)); +elevationAnglesToSat = zeros(size(sampleTimes)); + +progressStartText = "Starting ..."; +f = waitbar(0, progressStartText, "Name", "Task 1"); + +progressTotal = numel(sampleTimes); +progressQueue = parallel.pool.DataQueue; + +parfor_waitbar(f, progressTotal, progressStartText, true); +progressListener = afterEach(progressQueue, @(info) parfor_waitbar(f, progressTotal, info, false)); + +parfor idx = 1:numel(sampleTimes) + res = perform_access_analysis(satSim, sampleTimes(idx)); + distancesToSat(idx) = res.distanceToSatellite; + elevationAnglesToSat(idx) = res.elevationAngleToSatellite; + send(progressQueue, sprintf("Index %d", idx)); +end + +delete(progressListener); +if ishandle(f) + close(f); +end + +figure; +plot(sampleTimes, distancesToSat); +grid on; +xlabel("Time [s]"); +ylabel("Distance to Satellite [m]"); +title("Distance over time plot of satellite pass"); + +figure; +plot(sampleTimes, elevationAnglesToSat); +grid on; +xlabel("Time [s]"); +ylabel("Elevation angle to satellite [degree]"); +title("Elevation angle over time plot of satellite pass"); + +accessAnalysisResults.times = sampleTimes; +accessAnalysisResults.distances = distancesToSat; +accessAnalysisResults.elevationAngles = elevationAnglesToSat; \ No newline at end of file diff --git a/lab_5/task2.m b/lab_5/task2.m new file mode 100644 index 0000000..eea3de7 --- /dev/null +++ b/lab_5/task2.m @@ -0,0 +1,64 @@ +profile = "C"; +bandwidth = 20; +include_free_space_pathloss = true; +delay_model = "None"; +carrier_frequency = 2.4e9; +mobile_speed = 0; +mobile_altitude = groundStation.Altitude; +satellite_altitude = satConfig.satAltitude; +random_stream = "mt19937ar with seed"; +seed = 666; +tx_antenna_count = 1; +rx_antenna_count = 1; +apep = 1000; +txPower = 20; +wifiSensitivity = -85; +reliableRxPower = -75; +amplifierGain = 30; + +[cfgEHT, simParameters, channel, chInfo, maxChDelay, satelliteDopplerShift, chanBW] = create_baseline_configuration(profile, bandwidth, include_free_space_pathloss, delay_model, carrier_frequency, mobile_speed, mobile_altitude, satellite_altitude, random_stream, seed, tx_antenna_count, rx_antenna_count, apep); +cfgEHT.User{1}.MCS = 3; +simParameters.TransmitPower = txPower; + +[~, numPacketsSimulated, ~, rxPower, pathLoss] = simulateTransmission(accessAnalysisResults, cfgEHT, channel, maxChDelay, [], inf, simParameters); + +packetIdx = 1:numPacketsSimulated; +packetTimes = accessAnalysisResults.times(packetIdx); +rxPower = rxPower(packetIdx); +pathLoss = pathLoss(packetIdx); + +figure; +plot(packetTimes, rxPower); +grid on; +xlabel("Time"); +ylabel("Power [dBm]"); +title("Receive power over time plot of satellite pass"); +yline(wifiSensitivity, "--", "receive sensitivity -85 dBm"); +yline(reliableRxPower, "--", "reliable reception -75 dBm"); + +lambda = physconst("lightspeed") / carrier_frequency; +maxPathLossDetectable = txPower - wifiSensitivity; +maxDistanceDetectable = lambda / (4 * pi) * 10^(maxPathLossDetectable / 20); + +maxPathLossReliable = txPower - reliableRxPower; +maxDistanceReliable = lambda / (4 * pi) * 10^(maxPathLossReliable / 20); + +requiredTxPower = reliableRxPower + pathLoss; + +figure; +plot(packetTimes, requiredTxPower); +grid on; +xlabel("Time"); +ylabel("Power [dBm]"); +title("Required transmit power for -75 dBm receive power"); + +amplifiedTxPower = txPower + amplifierGain; +requiredAntennaGain = max(requiredTxPower - amplifiedTxPower, 0); + +validRequiredTxPower = requiredTxPower(isfinite(requiredTxPower)); +validRequiredAntennaGain = requiredAntennaGain(isfinite(requiredAntennaGain)); + +fprintf("Maximum detectable distance at %.0f dBm TX power: %.2f km\n", txPower, maxDistanceDetectable / 1e3); +fprintf("Maximum reliable distance at %.0f dBm TX power: %.2f km\n", txPower, maxDistanceReliable / 1e3); +fprintf("Required TX power for -75 dBm reception: %.2f dBm to %.2f dBm\n", min(validRequiredTxPower), max(validRequiredTxPower)); +fprintf("Required antenna gain with %.0f dB amplifier: %.2f dB to %.2f dB\n", amplifierGain, min(validRequiredAntennaGain), max(validRequiredAntennaGain)); diff --git a/tasks/05-leo-pass.pdf b/tasks/05-leo-pass.pdf new file mode 100644 index 0000000..75ddfd6 Binary files /dev/null and b/tasks/05-leo-pass.pdf differ