NR Positioning Using PRS

This example shows how to calculate the position of user equipment (UE) within a network of gNodeBs (gNBs) by using an NR positioning reference signal (PRS). The example uses the observed time difference of arrival (OTDOA) positioning approach to estimate the UE position in 2-D, relative to (0,0).

Introduction

The OTDOA positioning approach is one of the downlink-based UE positioning techniques. This technique uses the reference signal timing difference (RSTD) or time difference of arrival (TDOA) measurements to perform multilateration or trilateration by using the theory of hyperbolas. The RSTD is the relative timing difference between the neighboring gNB and the reference gNB (which might be the serving cell), as defined in TS 38.215 Section 5.1.29 [1]. For the OTDOA technique, the PRS is transmitted by a set of neighboring gNBs along with the serving gNB. In this case, one gNB acts as a reference cell for the RSTD measurements. To avoid a synchronization error in RSTD measurements, all of the gNBs transmit the PRS at the same time (which means synchronized transmission from all gNBs). The time of arrival (TOA) of the PRS signals from different gNBs are estimated at the UE side. RSTD values are computed for a set of neighboring gNB and reference gNB pairs. An RSTD value between neighboring ${\mathrm{gNB}}_{\mathit{j}}$ and reference ${\mathrm{gNB}}_{\mathit{i}}$ is defined as, ${\mathrm{RSTD}}_{\mathit{j},\mathit{i}}={\mathrm{TOA}}_{\mathit{j}}-{\mathrm{TOA}}_{\mathit{i}}$.

Each RSTD value that is obtained for a pair of gNBs corresponds to a hyperbola equation on which the UE is assumed to be located, with foci located at these gNBs. A set of hyperbola equations are deduced from a set of RSTD values, and the equations are solved to estimate the UE position. This process is called multilateration. The multilateration process involving one reference gNB and two neighboring gNBs is called trilateration (that is, the case of minimum requirement for UE positioning in 2-D).

This figure shows the positioning scenario with one reference gNB (which is the serving gNB in this case) and two neighboring gNBs. The PRS signals that are transmitted by all of the gNBs are received at the UE side at different timing instances (${\mathrm{TOA}}_{0}$, ${\mathrm{TOA}}_{1}$, and ${\mathrm{TOA}}_{2}$) based on the geographic distance between the UE and gNBs. From this network of three gNBs, two RSTD values are obtained from the TOA values. Each RSTD value corresponds to a hyperbola equation, on which the UE is assumed to be located. You can solve these two hyperbola equations to get UE position in 2-D.

This example shows how to create several gNB transmissions combined with different delays and received powers to model the reception of the waveforms from all of the gNBs by one UE in urban macro ('UMa') path loss scenario. The UE performs correlation with the PRS to establish the delay from each gNB and subsequently the delay difference between the neighboring gNBs and reference gNB. These delay differences (RSTD values) are used to compute the hyperbolas of constant delay difference and are plotted relative to known gNB positions. The intersection of these hyperbolas is computed to estimate the position of the UE. This example focuses on the 2-D positioning estimate of a UE by assuming that all gNBs are synchronized. This example gives a better estimate for UE position when the UE is located near the center of gNBs network. This example assumes that all gNBs are operated at the same carrier frequency (which corresponds to a single positioning frequency layer). This example also assumes that all of the carriers have the same subcarrier spacing, cyclic prefix, and frequency allocation within the common resource grid.

Simulation Parameters

Specify the number of frames and the carrier frequency.

```nFrames = 1; % Number of 10 ms frames fc = 3e9; % Carrier frequency (Hz)```

Specify UE and gNB positions. This example configures the positions in 3-D for 'UMa' path loss scenario but the position estimation algorithm in this example computes the position in 2-D (which means the estimated height is considered as 0 meters).

```% Configure UE position in xyz-coordinate plane UEPos = [500 -20 2]; % In meters. As per TR 38.901 Table 7.4.1-1, % UE height for 'UMa' path loss scenario is >= 1.5 meters and <= 22.5 meters % Configure number of gNBs and locate them at random positions in % xyz-coordinate plane numgNBs = 5; rng('default'); % Set RNG state for repeatability gNBPos = getgNBPositions(numgNBs); % In meters. This local function assumes all gNB % heights to be 25m, as per TR 38.901 Table 7.4.1-1, % for 'UMa' path loss scenario % Plot UE and gNB positions plotgNBAndUEPositions(gNBPos,UEPos,1:numgNBs);```

Configuration Objects

Carrier Configuration

Configure the carriers for all of the gNBs with different physical layer cell identities.

```cellIds = randperm(1008,numgNBs) - 1; % Configure carrier properties carrier = repmat(nrCarrierConfig,1,numgNBs); for gNBIdx = 1:numgNBs carrier(gNBIdx).NCellID = cellIds(gNBIdx); end validateCarriers(carrier);```

PRS Configuration

Configure PRS parameters for all of the gNBs. Configure the PRS for different gNBs such that no overlap exists to avoid the problem of hearability. For more information on how to configure PRS, see NR Positioning Reference Signal. This example considers different slot offsets for the PRS from different gNBs to avoid the overlap among PRS signals. This overlapping can be avoided in multiple ways by configuring the time-frequency aspects of PRS signals in an appropriate way (like choosing the nonoverlapped symbol allocation or frequency allocation, or by using muting pattern configuration parameters).

```% Slot offsets of different PRS signals prsSlotOffsets = 0:2:(2*numgNBs - 1); prsIDs = randperm(4096,numgNBs) - 1; % Configure PRS properties prs = nrPRSConfig; prs.PRSResourceSetPeriod = [10 0]; prs.PRSResourceOffset = 0; prs.PRSResourceRepetition = 1; prs.PRSResourceTimeGap = 1; prs.MutingPattern1 = []; prs.MutingPattern2 = []; prs.NumRB = 52; prs.RBOffset = 0; prs.CombSize = 12; prs.NumPRSSymbols = 12; prs.SymbolStart = 0; prs = repmat(prs,1,numgNBs); for gNBIdx = 1:numgNBs prs(gNBIdx).PRSResourceOffset = prsSlotOffsets(gNBIdx); prs(gNBIdx).NPRSID = prsIDs(gNBIdx); end```

PDSCH Configuration

Configure the physical downlink shared channel (PDSCH) parameters for all of the gNBs for the data transmission. This example assumes that the data is transmitted from all of the gNBs that have a single transmission layer.

```pdsch = nrPDSCHConfig; pdsch.PRBSet = 0:51; pdsch.SymbolAllocation = [0 14]; pdsch.DMRS.NumCDMGroupsWithoutData = 1; pdsch = repmat(pdsch,1,numgNBs); validateNumLayers(pdsch);```

Path Loss Configuration

Create a path loss configuration object for the `'Uma'` scenario.

```plCfg = nrPathLossConfig; plCfg.Scenario = 'Uma';```

Generate PRS and PDSCH Resources

Generate PRS and PDSCH resources corresponding to all of the gNBs. To improve the hearability of PRS signals, transmit PDSCH resources in the slots in which the PRS is not transmitted by any of the gNBs. This example generates PRS and PDSCH resources with unit power (0 dB).

```totSlots = nFrames*carrier(1).SlotsPerFrame; prsGrid = cell(1,numgNBs); dataGrid = cell(1,numgNBs); for slotIdx = 0:totSlots-1 [carrier(:).NSlot] = deal(slotIdx); [prsSym,prsInd] = deal(cell(1,numgNBs)); for gNBIdx = 1:numgNBs % Create an empty resource grid spanning one slot in time domain slotGrid = nrResourceGrid(carrier(gNBIdx),1); % Generate PRS symbols and indices prsSym{gNBIdx} = nrPRS(carrier(gNBIdx),prs(gNBIdx)); prsInd{gNBIdx} = nrPRSIndices(carrier(gNBIdx),prs(gNBIdx)); % Map PRS resources to slot grid slotGrid(prsInd{gNBIdx}) = prsSym{gNBIdx}; prsGrid{gNBIdx} = [prsGrid{gNBIdx} slotGrid]; end % Transmit data in slots in which the PRS is not transmitted by any of % the gNBs (to control the hearability problem) for gNBIdx = 1:numgNBs dataSlotGrid = nrResourceGrid(carrier(gNBIdx),1); if all(cellfun(@isempty,prsInd)) % Generate PDSCH indices [pdschInd,pdschInfo] = nrPDSCHIndices(carrier(gNBIdx),pdsch(gNBIdx)); % Generate random data bits for transmission data = randi([0 1],pdschInfo.G,1); % Generate PDSCH symbols pdschSym = nrPDSCH(carrier(gNBIdx),pdsch(gNBIdx),data); % Generate demodulation reference signal (DM-RS) indices and symbols dmrsInd = nrPDSCHDMRSIndices(carrier(gNBIdx),pdsch(gNBIdx)); dmrsSym = nrPDSCHDMRS(carrier(gNBIdx),pdsch(gNBIdx)); % Map PDSCH and its associated DM-RS to slot grid dataSlotGrid(pdschInd) = pdschSym; dataSlotGrid(dmrsInd) = dmrsSym; end dataGrid{gNBIdx} = [dataGrid{gNBIdx} dataSlotGrid]; end end % Plot carrier grid plotGrid(prsGrid,dataGrid);```

Perform OFDM Modulation

Perform orthogonal frequency-division multiplexing (OFDM) modulation of the signals at each gNB.

```% Perform OFDM modulation of PRS and data signal at each gNB txWaveform = cell(1,numgNBs); for waveIdx = 1:numgNBs carrier(waveIdx).NSlot = 0; txWaveform{waveIdx} = nrOFDMModulate(carrier(waveIdx),prsGrid{waveIdx} + dataGrid{waveIdx}); end % Compute OFDM information using first carrier, assuming all carriers are % at same sampling rate ofdmInfo = nrOFDMInfo(carrier(1));```

Add Signal Delays and Apply Path Loss

Calculate the time delays from each gNB to UE by using the known gNB and UE positions. This calculation uses the distance between the UE and gNB, `radius`, and the speed of propagation (speed of light). Calculate the sample delay by using the sampling rate, `info.SampleRate`, and store it in `sampleDelay`. This example only applies the integer sample delay by rounding the actual delays to their nearest integers. The example uses these variables to model the environment between gNBs and the UE. The information about these variables is not provided to the UE.

```speedOfLight = physconst('LightSpeed'); % Speed of light in m/s sampleDelay = zeros(1,numgNBs); radius = cell(1,numgNBs); for gNBIdx = 1:numgNBs radius{gNBIdx} = rangeangle(gNBPos{gNBIdx}',UEPos'); delay = radius{gNBIdx}/speedOfLight; % Delay in seconds sampleDelay(gNBIdx) = round(delay*ofdmInfo.SampleRate); % Delay in samples end```

Model the received signal at the UE by delaying each gNB transmission according to the values in `sampleDelay` and by attenuating the received signal from each gNB. Ensure that all of the waveforms are of the same length by padding the received waveform from each gNB with relevant number of zeros.

```rxWaveform = zeros(length(txWaveform{1}) + max(sampleDelay),1); rx = cell(1,numgNBs); for gNBIdx = 1:numgNBs % Calculate path loss for each gNB and UE pair losFlag = true; % Assuming the line of sight (LOS) flag as true, as we are only considering the LOS path delays in this example PLdB = nrPathLoss(plCfg,fc,losFlag,gNBPos{gNBIdx}(:),UEPos(:)); if PLdB < 0 || isnan(PLdB) || isinf(PLdB) error('nr5g:invalidPL',"Computed path loss (" + num2str(PLdB) + ... ") is invalid. Try changing the UE or gNB positions, or path loss configuration."); end PL = 10^(PLdB/10); % Add delay, pad, and attenuate rx{gNBIdx} = [zeros(sampleDelay(gNBIdx),1); txWaveform{gNBIdx}; ... zeros(max(sampleDelay)-sampleDelay(gNBIdx),1)]/sqrt(PL); % Sum waveforms from all gNBs rxWaveform = rxWaveform + rx{gNBIdx}; end```

TOA Estimation

Configure the number of cells or gNBs to be detected for the multilateration process. This example does not consider the transmission of primary and secondary synchronization signals for the purpose of cell search. Instead, the UE performs the correlation on the incoming signal with reference PRS generated for each gNB, and `cellsToBeDetected` number of best cells are selected based on correlation outcome. Compute TOAs of the signals from each gNB by using `nrTimingEstimate` function. To estimate the position of UE, compute RSTD values by using the absolute TOAs corresponding to the best cells.

```cellsToBeDetected = min(3,numgNBs); if cellsToBeDetected < 3 || cellsToBeDetected > numgNBs error('nr5g:InvalidNumDetgNBs',['The number of detected gNBs (' num2str(cellsToBeDetected) ... ') must be greater than or equal to 3 and less than or equal to the total number of gNBs (' num2str(numgNBs) ').']); end corr = cell(1,numgNBs); delayEst = zeros(1,numgNBs); maxCorr = zeros(1,numgNBs); for gNBIdx = 1:numgNBs [~,mag] = nrTimingEstimate(carrier(gNBIdx),rxWaveform,prsGrid{gNBIdx}); % Extract correlation data samples spanning about 1/14 ms for normal % cyclic prefix and about 1/12 ms for extended cyclic prefix (this % truncation is to ignore noisy side lobe peaks in correlation outcome) corr{gNBIdx} = mag(1:(ofdmInfo.Nfft*carrier(1).SubcarrierSpacing/15)); % Delay estimate is at point of maximum correlation maxCorr(gNBIdx) = max(corr{gNBIdx}); delayEst(gNBIdx) = find(corr{gNBIdx} == maxCorr(gNBIdx),1)-1; end % Get detected gNB numbers based on correlation outcome [~,detectedgNBs] = sort(maxCorr,'descend'); detectedgNBs = detectedgNBs(1:cellsToBeDetected); % Plot PRS correlation results plotPRSCorr(carrier,corr,ofdmInfo.SampleRate);```

UE Position Estimation Using OTDOA Technique

Calculate TDOA or RSTD values for all pairs of gNBs using the estimated TOA values by using the `getRSTDValues` function. Each RSTD value corresponding to a pair of gNBs can result from the UE being located at any position on a hyperbola with foci located at these gNBs.

```% Compute RSTD values for multilateration or trilateration rstdVals = getRSTDValues(delayEst,ofdmInfo.SampleRate);```

Consider RSTD values corresponding to the detected gNBs for the estimation of the UE position. This example assumes the first detected gNB as the reference gNB and the others as neighboring gNBs. Calculate a set of hyperbola equations by using RSTD values corresponding to the pairs of reference gNB and neighboring gNBs. Plot the hyperbola equations, where the intersection of these curves represents the UE position. If you want to know which gNBs correspond to a hyperbola curve, you can know that information from the data tip in the figure.

```% Plot gNB and UE positions txCellIDs = [carrier(:).NCellID]; cellIdx = 1; curveX = {}; curveY = {}; for jj = detectedgNBs(1) % Assuming first detected gNB as reference gNB for ii = detectedgNBs(2:end) rstd = rstdVals(ii,jj)*speedOfLight; % Delay distance % Establish gNBs for which delay distance is applicable by % examining detected cell identities txi = find(txCellIDs == carrier(ii).NCellID); txj = find(txCellIDs == carrier(jj).NCellID); if (~isempty(txi) && ~isempty(txj)) % Get x- and y- coordinates of hyperbola curve and store them [x,y] = getRSTDCurve(gNBPos{txi},gNBPos{txj},rstd); if isreal(x) && isreal(y) curveX{1,cellIdx} = x; %#ok<*SAGROW> curveY{1,cellIdx} = y; % Get the gNB numbers corresponding to the current hyperbola % curve gNBNums{cellIdx} = [jj ii]; cellIdx = cellIdx + 1; else warning("Due to the error in timing estimations, " + ... "hyperbola curve cannot be deduced for the combination of gNB" + jj + " and gNB" + ii+". So, ignoring this gNB pair for the estimation of UE position."); end end end end numHyperbolaCurves = numel(curveX); if numHyperbolaCurves < 2 error('nr5g:InsufficientCurves',"The number of hyperbola curves ("+numHyperbolaCurves+") is not sufficient for the estimation of UE position in 2-D. Try with more number of gNBs."); end ```

Use `getEstimatedUEPosition` function to estimate the UE position from the hyperbola curves. In general, two hyperbolas can intersect at atmost two points. If there are two points of intersection, the algorithm in this example returns the average of two points of intersection for UE position estimation.

```estimatedPos = getEstimatedUEPosition(curveX,curveY); % Compute positioning estimation error EstimationErr = norm(UEPos-estimatedPos); % in meters disp(['Estimated UE Position : [' num2str(estimatedPos(1)) ' ' num2str(estimatedPos(2)) ']' 13 ... 'UE Position Estimation Error: ' num2str(EstimationErr) ' meters']);```
```Estimated UE Position : [501.0045 -24.028] UE Position Estimation Error: 4.608 meters ```
```% Plot UE, gNB positions, and hyperbola curves gNBsToPlot = unique([gNBNums{:}],'stable'); plotPositionsAndHyperbolaCurves(gNBPos,UEPos,gNBsToPlot,curveX,curveY,gNBNums,estimatedPos);```

Summary and Further Exploration

This example shows OTDOA-based UE positioning estimation in 2-D.

You can vary the number of gNBs, locate the UE and gNBs at different positions.

You can also vary the time-frequency allocation parameters of PRS signals, vary the PDSCH configurations from each gNB, and see the variations in the estimated UE position.

References

1. 3GPP TS 38.215. "NR; Physical layer measurements (Release 16)." 3rd Generation Partnership Project; Technical Specification Group Radio Access Network.

2. 3GPP TR 38.901. “Study on channel model for frequencies from 0.5 to 100 GHz.” 3rd Generation Partnership Project; Technical Specification Group Radio Access Network.

Local Functions

This example uses these local functions to calculate gNB positions, validate carrier properties and the number of layers, estimate RSTD values, compute the UE position, and plot UE and gNB positions along with the generated hyperbolas.

```function gNBPos = getgNBPositions(numgNBs) % GNBPOS = getgNBPositions(NUMGNBS) returns a cell array of random % positions for gNBs, given the number of gNBs NUMGNBS. gNBPos = cell(1,numgNBs); for gNBIdx = 1:numgNBs % Position gNB randomly within gNBNum*2*pi/numgNBs radian sector phi = gNBIdx*2*pi/numgNBs + rand(1,1)*2*pi/(2*numgNBs) - 2*pi/(2*numgNBs); % Position gNB randomly between 4000 + (gNBNum*5000/numgNBs) and 5000 + % (gNBNum*5000/numgNBs) from UE r = randi([0,1000],1,1) + 4000 + (gNBIdx*5000/numgNBs); % Convert polar coordinates to Cartesian coordinates [x,y] = pol2cart(phi,r); gNBPos{gNBIdx} = [x,y,25]; % Considering all gNB heights to be 25m, as per TR 38.901 Table 7.4.1-1 for 'UMa' path loss scenario end end function plotgNBAndUEPositions(gNBPos,UEPos,gNBNums) % plotgNBAndUEPositions(GNBPOS,UEPOS,GNBNUMS) plots gNB and UE positions. numgNBs = numel(gNBNums); colors = getColors(numel(gNBPos)); figure; hold on; legendstr = cell(1,numgNBs); % Plot position of each gNB for gNBIdx = 1:numgNBs plot3(gNBPos{gNBNums(gNBIdx)}(1),gNBPos{gNBNums(gNBIdx)}(2),gNBPos{gNBNums(gNBIdx)}(3), ... 'Color',colors(gNBNums(gNBIdx),:),'Marker','^', ... 'MarkerSize',11,'LineWidth',2,'MarkerFaceColor',colors(gNBNums(gNBIdx),:)); legendstr{gNBIdx} = sprintf('gNB%d',gNBIdx); end % Plot UE position scatter3(UEPos(1),UEPos(2),UEPos(3),80,'LineWidth',2,'MarkerFaceColor','w','MarkerEdgeColor','k','MarkerFaceAlpha',0.2);view(2) axis equal; legend([legendstr 'UE']); xlabel('X Position (meters)'); ylabel('Y Position (meters)'); end function validateCarriers(carrier) % validateCarriers(CARRIER) validates the carrier properties of all % gNBs. % Validate physical layer cell identities cellIDs = [carrier(:).NCellID]; if numel(cellIDs) ~= numel(unique(cellIDs)) error('nr5g:invalidNCellID','Physical layer cell identities of all of the carriers must be unique.'); end % Validate subcarrier spacings scsVals = [carrier(:).SubcarrierSpacing]; if ~isscalar(unique(scsVals)) error('nr5g:invalidSCS','Subcarrier spacing values of all of the carriers must be same.'); end % Validate cyclic prefix lengths cpVals = {carrier(:).CyclicPrefix}; if ~all(strcmpi(cpVals{1},cpVals)) error('nr5g:invalidCP','Cyclic prefix lengths of all of the carriers must be same.'); end % Validate NSizeGrid values nSizeGridVals = [carrier(:).NSizeGrid]; if ~isscalar(unique(nSizeGridVals)) error('nr5g:invalidNSizeGrid','NSizeGrid of all of the carriers must be same.'); end % Validate NStartGrid values nStartGridVals = [carrier(:).NStartGrid]; if ~isscalar(unique(nStartGridVals)) error('nr5g:invalidNStartGrid','NStartGrid of all of the carriers must be same.'); end end function validateNumLayers(pdsch) % validateNumLayers(PDSCH) validates the number of transmission layers, % given the array of physical downlink shared channel configuration % objects PDSCH. numLayers = [pdsch(:).NumLayers]; if ~all(numLayers == 1) error('nr5g:invalidNLayers',['The number of transmission layers ' ... 'configured for the data transmission must be 1.']); end end function plotPRSCorr(carrier,corr,sr) % plotPRSCorr(CARRIER,CORR,SR) plots PRS-based correlation for each gNB % CORR, given the array of carrier-specific configuration objects CARRIER % and the sampling rate SR. numgNBs = numel(corr); colors = getColors(numgNBs); figure; hold on; % Plot correlation for each gNodeB t = (0:length(corr{1}) - 1)/sr; legendstr = cell(1,2*numgNBs); for gNBIdx = 1:numgNBs plot(t,abs(corr{gNBIdx}), ... 'Color',colors(gNBIdx,:),'LineWidth',2); legendstr{gNBIdx} = sprintf('gNB%d (NCellID = %d)', ... gNBIdx,carrier(gNBIdx).NCellID); end % Plot correlation peaks for gNBIdx = 1:numgNBs c = abs(corr{gNBIdx}); j = find(c == max(c),1); plot(t(j),c(j),'Marker','o','MarkerSize',5, ... 'Color',colors(gNBIdx,:),'LineWidth',2); legendstr{numgNBs+gNBIdx} = ''; end legend(legendstr); xlabel('Time (seconds)'); ylabel('Absolute Value'); title('PRS Correlations for All gNBs'); end function rstd = getRSTDValues(toa,sr) % RSTD = getRSTDValues(TOA,SR) computes RSTD values, given the time of % arrivals TOA and the sampling rate SR. % Compute number of samples delay between arrivals from different gNBs rstd = zeros(length(toa)); for jj = 1:length(toa) for ii = 1:length(toa) rstd(ii,jj) = toa(ii) - toa(jj); end end % Get RSTD values in time rstd = rstd./sr; end function [x,y] = getRSTDCurve(gNB1,gNB2,rstd) % [X,Y] = getRSTDCurve(GNB1,GNB2,RSTD) returns the x- and y-coordinates % of a hyperbola equation corresponding to an RSTD value, given the pair % of gNB positions gNB1 and gNB2. % Calculate vector between two gNBs delta = gNB1 - gNB2; % Express distance vector in polar form [phi,r] = cart2pol(delta(1),delta(2)); rd = (r+rstd)/2; % Compute the hyperbola parameters a = (r/2)-rd; % Vertex c = r/2; % Focus b = sqrt(c^2-a^2); % Co-vertex hk = (gNB1 + gNB2)/2; mu = -2:1e-3:2; % Get x- and y- coordinates of hyperbola equation x = (a*cosh(mu)*cos(phi)-b*sinh(mu)*sin(phi)) + hk(1); y = (a*cosh(mu)*sin(phi)+b*sinh(mu)*cos(phi)) + hk(2); end function plotPositionsAndHyperbolaCurves(gNBPos,UEPos,detgNBNums,curveX,curveY,gNBNums,estPos) % plotPositionsAndHyperbolaCurves(GNBPOS,UEPOS,DETGNBNUMS,CURVEX,CURVEY,GNBNUMS,ESTPOS) % plots gNB, UE positions, and hyperbola curves by considering these % inputs: % GNBPOS - Positions of gNBs % UEPOS - Position of UE % DETGNBNUMS - Detected gNB numbers % CURVEX - Cell array having the x-coordinates of hyperbola curves % CURVEY - Cell array having the y-coordinates of hyperbola curves % GNBNUMS - Cell array of gNB numbers corresponding to all hyperbolas % ESTPOS - Estimated UE position plotgNBAndUEPositions(gNBPos,UEPos,detgNBNums); for curveIdx = 1:numel(curveX) curves(curveIdx) = plot(curveX{1,curveIdx},curveY{1,curveIdx}, ... '--','LineWidth',0.9,'Color','k','MarkerFaceColor','k'); end % Add estimated UE position to figure plot3(estPos(1),estPos(2),estPos(3),'+','MarkerSize',10, ... 'MarkerFaceColor','#D95319','MarkerEdgeColor','#D95319','LineWidth',2); % Add legend gNBLegends = cellstr(repmat("gNB",1,numel(detgNBNums)) + ... detgNBNums + [" (reference gNB)" repmat("",1,numel(detgNBNums)-1)]); legend([gNBLegends {'Actual UE Position'} repmat({''},1,numel(curveX)) {'Estimated UE Position'}]); for curveIdx = 1:numel(curves) % Create a data tip and add a row to the data tip to display % gNBs information for each hyperbola curve dt = datatip(curves(curveIdx)); row = dataTipTextRow("Hyperbola of gNB" + gNBNums{curveIdx}(1) + ... " and gNB" + gNBNums{curveIdx}(2),''); curves(curveIdx).DataTipTemplate.DataTipRows = [row; curves(curveIdx).DataTipTemplate.DataTipRows]; % Delete the data tip that is added above delete(dt); end end function colors = getColors(numgNBs) % COLORS = getColors(NUMGNBS) returns the RGB triplets COLORS for the % plotting purpose, given the number of gNBs NUMGNBS. % Form RGB triplets if numgNBs <= 10 colors = [0.8500, 0.3250, 0.0980; 0.9290, 0.6940, 0.1250; 0.4940, 0.1840, 0.5560; ... 0.4660, 0.6740, 0.1880; 0.3010, 0.7450, 0.9330; 0.6350, 0.0780, 0.1840; ... 0.9290, 0.9, 0.3; 0.9290, 0.5, 0.9; 0.6660, 0.3740, 0.2880;0, 0.4470, 0.7410]; else % Generate 30 more extra RGB triplets than required. It is to skip % the gray and white shades as these shades are reserved for the % data and no transmissions in carrier grid plot in this example. % With this, you can get unique colors for up to 1000 gNBs. colors = colorcube(numgNBs+30); colors = colors(1:numgNBs,:); end end function plotGrid(prsGrid,dataGrid) % plotGrid(PRSGRID,DATAGRID) plots the carrier grid from all gNBs. numgNBs = numel(prsGrid); figure() mymap = [1 1 1; ... % White color for background 0.8 0.8 0.8; ... % Gray color for PDSCH data from all gNBs getColors(numgNBs)]; chpval = 3:numgNBs+2; gridWithPRS = zeros(size(prsGrid{1})); gridWithData = zeros(size(dataGrid{1})); names = cell(1,numgNBs); for gNBIdx = 1:numgNBs gridWithPRS = gridWithPRS + abs(prsGrid{gNBIdx})*chpval(gNBIdx); gridWithData = gridWithData + abs(dataGrid{gNBIdx}); names{gNBIdx} = ['PRS from gNB' num2str(gNBIdx)]; end names = [{'Data from all gNBs'} names]; % Replace all zero values of gridWithPRS with proper scaling (value of 1) % for white background gridWithPRS(gridWithPRS == 0) = 1; gridWithData(gridWithData ~=0) = 1; % Plot grid image(gridWithPRS+gridWithData); % In the resultant grid, 1 represents % the white background, 2 (constitute % of 1s from gridWithPRS and % gridWithData) represents PDSCH, and % so on % Apply colormap colormap(mymap); axis xy; % Generate lines L = line(ones(numgNBs+1),ones(numgNBs+1),'LineWidth',8); % Index colormap and associate selected colors with lines set(L,{'color'},mat2cell(mymap(2:numgNBs+2,:),ones(1,numgNBs+1),3)); % Set the colors % Create legend legend(names{:}); % Add title title('Carrier Grid Containing PRS and PDSCH from Multiple gNBs'); % Add labels to x-axis and y-axis xlabel('OFDM Symbols'); ylabel('Subcarriers'); end function estPos = getEstimatedUEPosition(xCell,yCell) % ESTPOS = getEstimatedUEPosition(XCELL,YCELL) returns the x-, y-, and z- % coordinates of the UE position, given the hyperbolas XCELL and YCELL. % Note that the z- coordinate is not estimated using this function but it % is returned as 0 meters. % Steps involved in this computation: % 1. Find closest points between hyperbolic surfaces % 2. Linearize surfaces around the points closest to intersection to % interpolate actual intersection location numCurves = numel(xCell); % Make all vectors have equal length maxLen = max(cellfun(@length,xCell)); for curveIdx = 1:numCurves xCell{curveIdx} = [xCell{curveIdx} inf(1,maxLen-length(xCell{curveIdx}))]; yCell{curveIdx} = [yCell{curveIdx} inf(1,maxLen-length(yCell{curveIdx}))]; end tempIdx = 1; for idx1 = 1:numCurves-1 for idx2 = (idx1+1):numCurves [firstCurve,secondCurve] = findMinDistanceElements(xCell{idx1},yCell{idx1},xCell{idx2},yCell{idx2}); for idx3 = 1:numel(firstCurve) [x1a,y1a,x1b,y1b] = deal(firstCurve{idx3}(1,1),firstCurve{idx3}(1,2), ... firstCurve{idx3}(2,1),firstCurve{idx3}(2,2)); [x2a,y2a,x2b,y2b] = deal(secondCurve{idx3}(1,1),secondCurve{idx3}(1,2), ... secondCurve{idx3}(2,1),secondCurve{idx3}(2,2)); a1 = (y1b-y1a)/(x1b-x1a); b1 = y1a - a1*x1a; a2 = (y2b-y2a)/(x2b-x2a); b2 = y2a - a2*x2a; xC(tempIdx) = (b2-b1)/(a1-a2); yC(tempIdx) = a1*xC(tempIdx) + b1; tempIdx = tempIdx+1; end end end estPos = [mean(xC) mean(yC) 0]; end function [firstCurvePoints,secondCurvePoints] = findMinDistanceElements(xA,yA,xB,yB) % [FIRSTCURVEPOINTS,SECONDCURVEPOINTS] = findMinDistanceElements(XA,YA,XB,YB) % returns the closest points between the given hyperbolic surfaces. distAB1 = zeros(numel(xA),numel(xB)); for idx1 = 1:numel(xA) distAB1(idx1,:) = sqrt((xB-xA(idx1)).^2 + (yB-yA(idx1)).^2); end [~,rows] = min(distAB1,[],'omitnan'); [~,col] = min(min(distAB1,[],'omitnan')); % Extract the indices of the points on two curves which are up to 5 % meters apart. This is to identify the multiple intersections between % two hyperbola curves. [allRows,allCols] = find(distAB1 <= distAB1(rows(col),col)+5); diffVec = abs(allRows - allRows(1)); index = find(diffVec > (min(diffVec) + max(diffVec))/2,1); allRows = [allRows(1);allRows(index)]; allCols = [allCols(1);allCols(index)]; for idx = 1:numel(allRows) firstCurveIndices = allRows(idx); secondCurveIndices = allCols(idx); x1a = xA(firstCurveIndices); y1a = yA(firstCurveIndices); x2a = xB(secondCurveIndices); y2a = yB(secondCurveIndices); % Use subsequent points to create line for linearization if firstCurveIndices == numel(rows) x1b = xA(firstCurveIndices-1); y1b = yA(firstCurveIndices-1); else x1b = xA(firstCurveIndices+1); y1b = yA(firstCurveIndices+1); end if secondCurveIndices == numel(rows) x2b = xB(secondCurveIndices-1); y2b = yB(secondCurveIndices-1); else x2b = xB(secondCurveIndices+1); y2b = yB(secondCurveIndices+1); end firstCurvePoints{idx} = [x1a y1a;x1b y1b]; %#ok<*AGROW> secondCurvePoints{idx} = [x2a y2a;x2b y2b]; end end```