5G Network Digital Twin in 3-D Urban Scenario
This example shows how to model a digital twin of a 5G network in a specific geographical location. The example performs a network or system-level simulation with ray-tracing-based channels.
Introduction
This example implements a 5G network simulation in a realistic urban environment using ray tracing for channel modeling. In this example, you create a digital twin of a network in a specific location and then evaluate its performance in that environment.
This example follows these steps:
The example models several base stations (gNBs) and user equipment (UE) units. You can configure the location of each gNB and UE. UEs remain static throughout the simulation. The ray tracing analysis results in propagation channel models for the network simulation.
After network simulation, the example calculates key performance indicators such as system throughput, spectral efficiency, and block error rate. It also provides UE-specific metrics including modulation and coding scheme (MCS) and signal to interference and noise ratio (SINR).
Specify System Parameters
Specify the main simulation parameters.
simulationTime = 0.025; % Seconds fc = 4.2e9; % Carrier frequency (Hz) PtxgNBdBm = 33; % gNB transmit power in dBm PtxUEdBm = 23; % UE transmit power in dBm NRB = 51; % Number of resource blocks SCS = 30; % Subcarrier spacing (kHz) duplexingMode = "TDD"; % "FDD" or "TDD" maxNumReflections = 2; % Maximum number of ray tracing reflections maxNumDiffractions = 1; % Maximum number of ray tracing diffractions gNBArrayRectangularDim = [2 2]; % gNB rectangular array dimensions (rows, columns). Valid values: % [1 1], [2 1], [2 2], [2 4], [4 4] ueArrayRectangularDim = [2 1]; % UE rectangular array dimensions (rows, columns) validategNBArraySize(gNBArrayRectangularDim); BW = validateNumerology(SCS,NRB);% BW: bandwidth in Hz
Specify the name of an OpenStreetMap file [1] with 3-D buildings information, and specify the geographic coordinates (latitude and longitude) of the gNBs and UEs.
fileName = "chicago.osm";The number of gNBs and UEs is specified by the size of the location arrays gNBLocations and ueLocations.
% gNB locations gNBLocations = [41.876853,-87.626942; % gNB1 41.877706,-87.630108; % gNB2 41.875433,-87.630880]; % gNB3 gNBAzAngle = [180 -90 70]; % gNB array normal azimuth angles (deg) relative to east gNBElAngle = [-6 -6 -6]; % gNB array normal elevation angles (deg), negative is downwards gNBAntennaHeight = [5 5 25]; % gNB antenna height (m) above mounting surface (ground or building roof) numgNBs = size(gNBLocations,1); gNBLabels = "gNB"+string(1:numgNBs); % UE locations ueLocations = [41.877306,-87.627886; % UE1 41.877090,-87.628397; % UE2 41.877448,-87.629231; % UE3 41.876888,-87.630645; % UE4 41.876920,-87.629331; % UE5 41.875769,-87.628913; % UE6 41.876582,-87.630087]; % UE7 ueAntennaHeight = 1.5; % m numUEs = size(ueLocations,1); ueLabels = "UE"+string(1:numUEs); % Random number generator default settings for result reproducibility rng("default");
Display the gNB and UE locations.
plotNodeLocations(fileName,"gNB and UE Locations", ... gNBLocations,ueLocations,gNBLabels,ueLabels);
Perform Ray Tracing Analysis
Perform ray tracing analysis, resulting in a set of rays the example uses for channel modeling during network simulation.
Create gNB Sites
Create the gNB sites for ray tracing analysis.
Precalculate the gNB transmit power in watts, the wavelength, and the reference ellipsoid for converting geographical coordinates (latitude and longitude) to cartesian coordinates used by the network simulation nodes.
% gNB transmit power in watts PtxW = db2pow(PtxgNBdBm-30); % Wavelength for array dimensioning lambda = freq2wavelen(fc); % Reference ellipsoid object for lat/lon to cartesian conversion wgs84 = wgs84Ellipsoid;
Create the gNB site objects for ray tracing analysis. Each gNB has its own array with the broadside pointing in the direction specified in gNBAzAngle and gNBElAngle. The gNBsTable table contains name, indexing, and cartesian coordinates information associated with each gNB. You can use this information later to create the nodes for network simulation.
% Create gNB site objects (txsite objects) gNBsTable = table(Size=[numgNBs 4], ... VariableTypes=["string","int8","int8","double"], ... VariableNames=["gNBLabel","gnbArrayIdx","gnbID","xyz"]); gNBsTable.xyz = zeros(numgNBs,3); % Make each xyz row a 3-element vector: [x y z] gnbSites = txsite.empty(); % Calculate bandwidth and uplink and downlink carrier frequencies [dlFc,ulFc] = fddCarrierFreq(duplexingMode,fc,BW); for n = 1:numgNBs % Array object txArray = phased.NRRectangularPanelArray(Size=[gNBArrayRectangularDim 1 1], ... Spacing=[0.5*lambda*[1 1] 1 1], ... ElementSet={phased.NRAntennaElement(PolarizationAngle=45), ... phased.NRAntennaElement(PolarizationAngle=-45)}); % Create gNB txsites gnbSites(n)= txsite(Name=(gNBLabels(n)), ... Latitude=gNBLocations(n,1), ... Longitude=gNBLocations(n,2), ... TransmitterPower=PtxW, ... TransmitterFrequency=dlFc, ... Antenna=txArray, ... AntennaHeight=gNBAntennaHeight(n), ... AntennaAngle=[gNBAzAngle(n) gNBElAngle(n)]'); % Equivalent cartesian coordinates: gNB of interest is origin [x,y,z] = latLonToCartesian(gNBLocations(n,:),gnbSites(n).elevation+gnbSites(n).AntennaHeight,gNBLocations(1,:),wgs84); gNBsTable(n,:) = {gNBLabels(n),n,0,[x,y,z]}; end % Assume all gNBs have the same number of antennas numgNBAntennas = txArray.getNumElements;
Create UE Sites
Create the UE site objects for ray tracing analysis. The uesTable table contains name, indexing, and cartesian coordinates information associated with each UE. You can use this information later to create the nodes for network simulation.
% Create UE dipole antenna elements pol1AntElement = createUEAntennaElement(0,fc); % Antenna element with 1st polarization angle pol2AntElement = createUEAntennaElement(90,fc); % Antenna element with 2nd polarization angle % Table with UE information uesTable = table(Size=[numUEs 5], ... VariableTypes=["string","int8","int8","double","double"], ... VariableNames=["ueLabel","ueArrayIdx","ueID","latLon","xyz"]); uesTable.latLon = zeros(numUEs,2); % Make each latLon row a 2-element vector: [lat lon] uesTable.xyz = zeros(numUEs,3); % Make each xyz row a 3-element vector: [x y z] % Generate UE sites ueSites = rxsite.empty(); for n = 1:numUEs rxArray = phased.NRRectangularPanelArray(Size=[ueArrayRectangularDim 1 1], ... Spacing=[0.5*lambda*[1 1] 1 1], ... ElementSet={clone(pol1AntElement),clone(pol2AntElement)}); % Create UE rxsites ueSites(n) = rxsite(Name=ueLabels(n), ... Latitude=ueLocations(n,1), ... Longitude=ueLocations(n,2), ... Antenna=rxArray, ... AntennaHeight=ueAntennaHeight); % Equivalent cartesian coordinates: gNB of interest is origin [x,y,z] = latLonToCartesian(ueLocations(n,:),ueSites(n).elevation+ueSites(n).AntennaHeight,gNBLocations(1,:),wgs84); uesTable(n,:) = {ueLabels(n),n,0,[ueLocations(n,:)],[x,y,z]}; end % Assume all UEs have the same number of antennas numUEAntennas = rxArray.getNumElements;
Open Site Viewer to view the map and buildings, and display the gNB and UE locations.
% Open Site viewer and display tx and rx sites if exist('viewer','var') && isvalid(viewer) % viewer handle exists and viewer window is open clearMap(viewer); else viewer = siteviewer(Basemap="openstreetmap",Buildings=fileName); end show(ueSites) show(gnbSites)
Perform Ray Tracing Analysis
Perform ray tracing analysis between each gNB and all UEs. Create a ray tracing propagation model with the specified maximum number of reflections and diffractions. You can specify other parameters associated with ray tracing in this propagation model object.
% Ray tracing propagation model pm = propagationModel("raytracing",MaxNumReflections=maxNumReflections,MaxNumDiffractions=maxNumDiffractions,MaxRelativePathLoss=Inf);
The gnbSites and ueSites are arrays of txsite and rxsite objects respectively. You can use them to perform the downlink ray tracing analysis. To perform the uplink ray tracing analysis, create txsite and rxsite objects for UEs and gNBs respectively. You can also use these objects to create the uplink ray tracing channel models.
% Create sites for UL modeling
[txueUL,rxgNBUL] = createULSites(gnbSites,ueSites,db2pow(PtxUEdBm-30),ulFc);Perform the ray tracing analysis between each gNB and each UE in the uplink (UL) and downlink (DL) direction.
For TDD, the example assumes channel reciprocity. In that case, it skips the uplink ray tracing calculation. Instead, it uses the information in the downlink rays to calculate the uplink rays.
raysDL = cell(1,numgNBs); raysUL = cell(1,numgNBs); for n = 1:numgNBs % DL ray tracing analysis raysDL{n} = raytrace(gnbSites(n),ueSites,pm); if duplexingMode == "TDD" % Assume channel reciprocity, no need to ray trace in the uplink, raysUL{n} = dl2ulRays(raysDL{n}); else % FDD % Uplink ray tracing analysis for m=1:numel(txueUL) ulrays = raytrace(txueUL(m),rxgNBUL(n),pm); raysUL{n}{m} = ulrays{1}; end end end
Simulate Network
Now set up and perform the network simulation using the results of the ray tracing analysis.
Initialize the network simulator.
networkSimulator = wirelessNetworkSimulator.init;
Create gNB Nodes
In the previous sections you used txsite and rxsite objects to perform ray tracing analysis. The network simulation uses nrGNB and nrUE objects to represent the different nodes in the simulation.
Create nrGNB nodes for network simulation. Specify the position of each node using the cartesian coordinates that correspond to the geographical coordinates (latitude and longitude) stored in gNBsTable.
gnbsxyz = gNBsTable.xyz; gnbs = nrGNB(Position=gnbsxyz, ... NumTransmitAntennas=numgNBAntennas(1), ... NumReceiveAntennas=numgNBAntennas(1), ... CarrierFrequency=fc, ... ChannelBandwidth=BW, ... NumResourceBlocks=NRB, ... SubcarrierSpacing=SCS, ... DuplexMode=duplexingMode, ... TransmitPower=PtxgNBdBm, ... NoiseFigure=6, ... ReceiveGain=6, ... SRSPeriodicityUE=10, ... PHYModel="abstract-phy", ... Name=gNBsTable.gNBLabel'); % gNB node IDs to table gNBsTable.gnbID = [gnbs.ID]';
Configure the scheduler. Specify the scheduling algorithm, the downlink channel state information measurement signal, and the link adaptation configuration parameters.
% Link adaptation structure laConfig = struct(StepUp=0.9,StepDown=0.1,InitialOffset=1,ResetOffsetOnCSIReport=false); % Scheduler configuration configureScheduler(gnbs,Scheduler="RoundRobin",CSIMeasurementSignalDL="CSI-RS", ... LinkAdaptationConfigDL=laConfig,LinkAdaptationConfigUL=laConfig);
Create UE Nodes
Create nrUE nodes per UE for network simulation. Specify the position of each node using the cartesian coordinates that correspond to the geographical coordinates (latitude and longitude) stored in uesTable.
% Create UE objects uesxyz = uesTable.xyz; ues = nrUE(Position=uesxyz, ... NumTransmitAntennas=numUEAntennas(1), ... NumReceiveAntennas=numUEAntennas(1), ... TransmitPower=PtxUEdBm, ... NoiseFigure=6, ... ReceiveGain=6, ... Name=uesTable.ueLabel'); % UE node IDs to table uesTable.ueID = [ues.ID]';
Set Up Network Simulator
Add the nrGNB and nrUE nodes to the network simulator.
addNodes(networkSimulator,gnbs) addNodes(networkSimulator,ues)
Connect UEs to gNBs
Connect the UE nodes to the gNB nodes and enable full-buffer traffic.
Generate a connection matrix where each row corresponds to a UE and each column corresponds to a gNB. An entry of 1 in the matrix indicates that the corresponding UE (row) is connected to the corresponding gNB (column). In this example, the system assigns each UE to the closest gNB.
[connectionMatrix,gNBsTable,uesTable] = generateConnectionMatrix(gnbSites,ueSites,gNBsTable,uesTable);
Connect gNBs to UEs using the connection matrix information. Specify full-buffer traffic for both the uplink and downlink direction.
for n=1:numgNBs gnbIdx = gNBsTable.gnbArrayIdx(n); targetUEsIdx = gNBsTable.AssociatedUEsIdx{n}; if ~isempty(targetUEsIdx) connectUE(gnbs(gnbIdx),ues(targetUEsIdx),FullBufferTraffic="on"); end associatedUEsStr = mat2str([ues(cell2mat(gNBsTable.AssociatedUEsIdx(n))).Name]); disp("UEs associated to gNB "+gnbs(n).ID+" ("""+gnbs(n).Name+"""): "+associatedUEsStr) end
UEs associated to gNB 1 ("gNB1"): ["UE1" "UE2"]
UEs associated to gNB 2 ("gNB2"): ["UE3" "UE4" "UE5" "UE7"]
UEs associated to gNB 3 ("gNB3"): "UE6"
Create Ray Tracing Channel Models
Use the ray tracing rays to create a ray tracing channel for each gNB-to-UE link.
Calculate OFDM-related information. The example uses this information later to obtain the signal sampling rate and related 5G numerology information.
% Calculate OFDM information waveformInfo = nrOFDMInfo(NRB,SCS); sr = waveformInfo.SampleRate; % Channel filtering flag based on PHY model used. Assume all gNBs use the % same PHY model channelFilteringFlag = (gnbs(1).PHYModel~="abstract-phy"); rtChannelRegistry = RayTracingChannelRegistry([gnbs.ID ues.ID],sr,channelFilteringFlag);
Create a ray tracing channel for each gNB-UE pair.
createRayTracingChannels(rtChannelRegistry,raysDL,raysUL, ... gnbSites,ueSites,txueUL,rxgNBUL, ... gnbs,ues);
Use a function handle to add the ray tracing channel model to the network simulator.
addChannelModel(networkSimulator,@rtChannelRegistry.applyChannel);
Run the Simulation
Create the metric visualizer and the simulation scheduling logger for the gNBs to store simulation metrics and run-time scheduling information.
[metricsVisualizer,simSchedulingLogger] = createVisualizerAndLogger(simulationTime,1:numgNBs,gnbs,ues);
Run the simulation for the specified simulation time.
run(networkSimulator,simulationTime)
Obtain the log tables for all the UEs.
logTables = cell(numgNBs,1); for n = 1:numgNBs uesOfInterest = ues(ismember([ues.ID],[gnbs(n).UENodeIDs])); % UEs associated to this gNB logTables{n} = getSimulationLogsTables(simSchedulingLogger{n},gnbs(n).DuplexMode,[uesOfInterest.RNTI]); end
Evaluate Simulation Results
You can inspect the simulation results using the generated log tables. This example shows how to obtain the downlink and uplink throughput, the average downlink and uplink SINR, and the MCS at the end of the simulation for all the gNBs and their associated UEs.
gNBsMetricsSummary = metricsSummary(gnbs,ues,metricsVisualizer,logTables,simulationTime); display(gNBsMetricsSummary)
gNBsMetricsSummary=7×10 table
gNB UE DL Thput (Mbps) Avg DL SINR (dB) DL MCS Max DL NumLayers UL Thput (Mbps) Avg UL SINR (dB) UL MCS Max UL NumLayers
______ _____ _______________ ________________ ______ ________________ _______________ ________________ ______ ________________
"gNB1" "UE1" 38.862 22.541 25 2 47.444 24.933 24 3
"gNB1" "UE2" 33.771 16.396 15 3 53.093 25.121 26 3
"gNB2" "UE3" 6.7888 13.478 21 2 20.445 27.539 26 2
"gNB2" "UE4" 0.04544 8.4092 22 2 13.318 16.502 18 4
"gNB2" "UE5" 2.4666 10.828 23 2 15.656 17.763 20 4
"gNB2" "UE7" 6.6589 15.415 20 2 23.525 21.474 21 4
"gNB3" "UE6" 23.34 13.289 18 2 57.873 23.767 24 2
Display downlink throughput per UE on the map.
figure; [~,loc] = ismember(gNBsMetricsSummary.UE, uesTable.ueLabel); plotNodeLocations(fileName,"Downlink Throughput (Mbps)",gNBLocations, ... uesTable.latLon(loc,:),gNBLabels,gNBsMetricsSummary.UE,gNBsMetricsSummary.("DL Thput (Mbps)"));
You can also display the peak and achieved throughput, spectral efficiency, and block error rate.
for m=1:numgNBs disp(gnbs(m).Name+"---------------------------"+newline) if ~isempty(metricsVisualizer{m}) displayPerformanceIndicators(metricsVisualizer{m}); end end
gNB1---------------------------
Peak UL throughput: 210.31 Mbps Achieved cell UL throughput: 100.54 Mbps Achieved UL throughput for each UE: [47.44 53.09] Peak UL spectral efficiency: 10.52 bits/s/Hz Achieved UL spectral efficiency for cell: 5.03 bits/s/Hz Block error rate for each UE in the UL direction: [0.222 0.111] Peak DL throughput: 290.09 Mbps Achieved cell DL throughput: 72.63 Mbps Achieved DL throughput for each UE: [38.86 33.77] Peak DL spectral efficiency: 14.50 bits/s/Hz Achieved DL spectral efficiency for cell: 3.63 bits/s/Hz Block error rate for each UE in the DL direction: [0.138 0.31]
gNB2---------------------------
Peak UL throughput: 210.31 Mbps Achieved cell UL throughput: 72.94 Mbps Achieved UL throughput for each UE: [20.44 13.32 15.66 23.53] Peak UL spectral efficiency: 10.52 bits/s/Hz Achieved UL spectral efficiency for cell: 3.65 bits/s/Hz Block error rate for each UE in the UL direction: [0 0.611 0.556 0.333] Peak DL throughput: 290.09 Mbps Achieved cell DL throughput: 15.96 Mbps Achieved DL throughput for each UE: [6.79 0.05 2.47 6.66] Peak DL spectral efficiency: 14.50 bits/s/Hz Achieved DL spectral efficiency for cell: 0.80 bits/s/Hz Block error rate for each UE in the DL direction: [0.621 0.897 0.793 0.586]
gNB3---------------------------
Peak UL throughput: 210.31 Mbps Achieved cell UL throughput: 57.87 Mbps Achieved UL throughput for each UE: [57.87] Peak UL spectral efficiency: 10.52 bits/s/Hz Achieved UL spectral efficiency for cell: 2.89 bits/s/Hz Block error rate for each UE in the UL direction: [0.222] Peak DL throughput: 290.09 Mbps Achieved cell DL throughput: 23.34 Mbps Achieved DL throughput for each UE: [23.34] Peak DL spectral efficiency: 14.50 bits/s/Hz Achieved DL spectral efficiency for cell: 1.17 bits/s/Hz Block error rate for each UE in the DL direction: [0.655]
Using the log tables generated, you can investigate how certain metrics evolve throughout the simulation. For example, you can inspect the MCS and SINR values in each slot for the UEs.
plotMCSandSINR(logTables,ues,gnbs,waveformInfo);
The simulation scenario is static. Therefore any variations in the results displayed are caused by the dynamic interference introduced by the scheduler, and the link adaptation algorithm. Depending on which UE is scheduled for downlink transmission, the interference level on other UEs can fluctuate. The downlink MCS decreases as the simulation progresses. This decrease is introduced by the link adaptation algorithm. This behavior happens because:
The initial channel state information (CSI) estimation uses reference signals without precoding (CSI-RS and SRS). The gNB chooses the channel rank and the MCS using the SINR measured during the CSI estimation.
Data packets use precoding. Therefore, the SINR experienced by UEs is different from the SINR observed during the CSI estimation.
This SINR mismatch can result in packet reception errors. The link adaptation algorithm adjusts the downlink MCS to compensate for the SINR mismatch, reduce the probability of packet reception failure, and increase the throughput.
References
[1] The OpenStreetMap file is downloaded from https://www.openstreetmap.org, which provides access to crowd-sourced map data all over the world. The data is licensed under the Open Data Commons Open Database License (ODbL), https://opendatacommons.org/licenses/odbl/.
Local Functions
function [x,y,z] = latLonToCartesian(latLon, totalAntennaHeight, latLonReference, refEllipsoid) % Convert lat/lon and height position into cartesian representation (x,y,z) % relative to a reference point. LATLON is a two-element vector with the % latitude and longitude values to convert, and TOTALANTENNAHEIGHT is the % total antenna height, including terrain, building, and antenna height. % LATLONREF is a two element vector with the reference latitude and % longitude. The conversion assumes a height of 0 meters for the reference % point. REFELLIPSOID is the reference ellipsoid object [x,y,z] = geodetic2enu(latLon(1),latLon(2),totalAntennaHeight,latLonReference(1),latLonReference(2),0,refEllipsoid); end function validategNBArraySize(arraySize) % Validate arraySize ([rows columns]) array size. % Valid sizes: % [1 1] % [2 1] % [2 2] % [2 4] % [4 4] % Valid sizes validSizes = [1 1; 2 1; 2 2; 2 4; 4 4]; % Check arraySize validity if ~ismember(arraySize(:).',validSizes,"rows") error("Invalid gNB array size ("+mat2str(arraySize)+"). Valid sizes are: "+mat2str(validSizes)) end end function BW = validateNumerology(SCS,NRB) % Validate subcarrier spacing and grid size values and return the bandwidth in Hz % Channel transmission bandwidth configuration table (SCS, BW and NRB) txBwConfig = nrDLCarrierConfig.FR1BandwidthTable; rowName = string(SCS+"kHz"); % SCS to string % Validate SCS if ~ismember(rowName, string(txBwConfig.Properties.RowNames)) error("Invalid SCS ("+SCS+"). Value must be a member of: "+ ... strjoin(erase(string(txBwConfig.Properties.RowNames),"kHz"))); end nrbRow = txBwConfig{rowName, :}; % Extract the row of NRB values validMask = ~isnan(nrbRow); % Valid entries are those that are not NaN matchMask = validMask & (nrbRow == NRB); % Find exact NRB matches among valid entries if any(matchMask) % Column names where there is a match bwVars = string(txBwConfig.Properties.VariableNames(matchMask)); % Convert "XXMHz" to numeric Hz BW = str2double(erase(bwVars, "MHz"))*1e6; % Remove the "MHz" suffix and convert to double else disp("Valid channel transmission bandwidth configurations: SCS, bandwidth and NRBs.") disp(txBwConfig) [~,idx]=min(abs(nrbRow-NRB)); error("Invalid combination of SCS ("+SCS+" kHz) and NRB ("+NRB+"). "+newline+ ... "Select a valid combination from the channel transmission bandwidth configuration table."+newline+ ... "The closest set of values in the table is SCS = "+SCS+"kHz, and NRB = "+nrbRow(idx)); end end function [dlFc,ulFc] = fddCarrierFreq(duplexingMode,fc,channelBandwidth) % Return the downlink and uplink channels center carrier frequency if duplexingMode=="TDD" dlFc = fc; ulFc = fc; else % DL/UL band offset is +/- guardBand*0.5 Hz of the carrier frequency and is % channelBandwidth Hz wide. guardBand = 140e6; dlFc = fc + guardBand/2 + channelBandwidth/2; ulFc = fc - guardBand/2 - channelBandwidth/2; end end function plotMCSandSINR(logTables,ues,gnbs,waveformInfo) % Plot the MCS and SINR per slot for all gNBs and UEs simulated. for m=1:numel(logTables) if ~isempty(logTables{m}) uesOfInterest = ues(ismember([ues.ID],[gnbs(m).UENodeIDs])); ueName = [uesOfInterest.Name]; figure; tl = tiledlayout(2,3); rntisLocal = unique(logTables{m}.grantLogsTable.RNTI); % UEs assigned to this gNB for n=1:numel(rntisLocal) [dlmcs,dlsinr,dlnumLayers,dlSlotNum] = getMCSSINR(logTables{m}.dlUELogs{n},waveformInfo.SlotsPerFrame); [ulmcs,ulsinr,ulnumLayers,ulSlotNum] = getMCSSINR(logTables{m}.ulUELogs{n},waveformInfo.SlotsPerFrame); nexttile(tl,1); plot(dlSlotNum,dlmcs,'.'); hold on; grid on; title("Downlink MCS"); ylabel("MCS") nexttile(tl,2); plot(dlSlotNum,dlsinr,'.'); hold on; grid on; title("Downlink SINR"); ylabel("SINR (dB)"); nexttile(tl,3); plot(dlSlotNum,dlnumLayers,'.'); hold on; grid on; title("Downlink Num Layers"); ylabel("Num Layers"); nexttile(tl,4); plot(ulSlotNum,ulmcs,'.'); hold on; grid on; title("Uplink MCS"); ylabel("MCS") nexttile(tl,5); plot(ulSlotNum,ulsinr,'.'); hold on; grid on; title("Uplink SINR"); ylabel("SINR (dB)") nexttile(tl,6); plot(ulSlotNum,ulnumLayers,'.'); hold on; grid on; title("Uplink Num Layers"); ylabel("Num Layers"); end legendtext = ueName; legend(legendtext) xlabel(tl,"Slot number"); title(tl,gnbs(m).Name+". MCS and SINR ("+gnbs(1).DuplexMode+")") % assume all gNBs have the same duplexing mode leg = legend(legendtext,Orientation="vertical",NumColumns=2);leg.Layout.Tile = 'east'; end end end function ueAntennaElement = createUEAntennaElement(polAngleDeg,fc) % Create a custom antenna element using a short dipole element as a % reference. POLANGLEDEG is the polarization angle in degrees relative to % the z axes of the reference short dipole element. FC is the frequency % used to calculate the antenna response. % Compute polarization direction vector polVector = [0; sind(polAngleDeg); cosd(polAngleDeg)]; % y–z plane rotation % Create reference short dipole antenna element shortDipole = phased.ShortDipoleAntennaElement(AxisDirection="Custom",CustomAxisDirection=polVector); % Define angles azAngles = -180:180; % Azimuth angles elAngles = -90:90; % Elevation angles % Extract H and V patterns in magnitude and phase patHmag = zeros(numel(elAngles),numel(azAngles)); patVmag = patHmag; patHphase = patHmag; patVphase = patHmag; for m = 1:numel(elAngles) pat = step(shortDipole,fc,[azAngles;elAngles(m)*ones(1,numel(azAngles))]); patHmag(m,:) = mag2db(abs(pat.H)); % Store horizontal pattern magnitude patVmag(m,:) = mag2db(abs(pat.V)); % Store vertical pattern magnitude patHphase(m,:) = rad2deg(angle(pat.H)); % Store horizontal pattern phase patVphase(m,:) = rad2deg(angle(pat.V)); % Store vertical pattern phase end % Initialize the custom antenna element with polarization applied ueAntennaElement = phased.CustomAntennaElement(PatternCoordinateSystem='az-el', ... AzimuthAngles=azAngles,ElevationAngles=elAngles, ... SpecifyPolarizationPattern=true, ... HorizontalMagnitudePattern=patHmag,HorizontalPhasePattern=patHphase, ... VerticalMagnitudePattern=patVmag,VerticalPhasePattern=patVphase); end function [metricsVisualizer,simSchedulingLogger] = createVisualizerAndLogger(simulationTime,gNBsOfInterestID,gnbs,ues) % Create metric visualizer and scheduling logger objects for the gNBs of % interest numFrameSimulation = simulationTime/10e-3; metricsVisualizer = cell(numel(gNBsOfInterestID),1); simSchedulingLogger = cell(numel(gNBsOfInterestID),1); for n = 1:numel(gNBsOfInterestID) uesOfInterest = ues(ismember([ues.ID],[gnbs(gNBsOfInterestID(n)).UENodeIDs])); % UEs associated to this gNB if ~isempty(uesOfInterest) metricsVisualizer{n} = helperNRMetricsVisualizer(gnbs(gNBsOfInterestID(n)),uesOfInterest, ... PlotSchedulerMetrics=false,PlotPhyMetrics=false); % Create an object to log scheduler traces simSchedulingLogger{n} = helperNRSchedulingLogger(numFrameSimulation,gnbs(gNBsOfInterestID(n)),uesOfInterest); else metricsVisualizer{n} = []; end end end function metricsSummary = metricsSummary(gnbs,ues,metricsVisualizer,logTables,simulationTime) % Compute and summarize key performance metrics for each UE associated with % gNBs of interest and aggregate results into a table. The summary includes % throughput, average SINR, and MCS for both downlink and uplink per UE. metricsSummary = []; numgNBs = numel(gnbs); for m=1:numgNBs % UEs associated to this gNB uesAssociatedTogNB = ues(ismember([ues.ID],[gnbs(m).UENodeIDs])); if ~isempty(uesAssociatedTogNB) rntis = zeros(numel(uesAssociatedTogNB),1); dlAvgSINR = zeros(numel(uesAssociatedTogNB),1); ulAvgSINR = zeros(numel(uesAssociatedTogNB),1); dlMCS = zeros(numel(uesAssociatedTogNB),1); ulMCS = zeros(numel(uesAssociatedTogNB),1); dlPeakNumLayers = zeros(numel(uesAssociatedTogNB),1); ulPeakNumLayers = zeros(numel(uesAssociatedTogNB),1); % Calculate average SINR and MCS slotOffset = 10; % start slot to obtain average metrics (avoid transient) % if simulation has less than 10 slots, use two thirds for % average metrics for n=1:numel(uesAssociatedTogNB) if slotOffset > numel(logTables{m}.dlUELogs{n}.SINR) slotOffset = floor(numel(logTables{m}.dlUELogs{n}.SINR)/3); end rntis(n) = uesAssociatedTogNB(n).RNTI; dlAvgSINR(n) = mean(logTables{m}.dlUELogs{n}.SINR(slotOffset:end)); ulAvgSINR(n) = mean(logTables{m}.ulUELogs{n}.SINR(slotOffset:end)); dlMCS(n) = logTables{m}.dlUELogs{n}.("MCS Index")(end); ulMCS(n) = logTables{m}.ulUELogs{n}.("MCS Index")(end); dlPeakNumLayers(n) = max(logTables{m}.dlUELogs{n}.NumLayers); ulPeakNumLayers(n) = max(logTables{m}.ulUELogs{n}.NumLayers); end [~, rntiIdx] = sort(rntis); dlAvgSINR = dlAvgSINR(rntiIdx); ulAvgSINR = ulAvgSINR(rntiIdx); dlMCS = round(dlMCS(rntiIdx)); ulMCS = round(ulMCS(rntiIdx)); dlPeakNumLayers = dlPeakNumLayers(rntiIdx); ulPeakNumLayers = ulPeakNumLayers(rntiIdx); % Calculate throughput simMetrics = getMetrics(metricsVisualizer{m}); dlTh = 1e-6*8*simMetrics.MACMetrics.("DL Throughput Bytes")/simulationTime; ulTh = 1e-6*8*simMetrics.MACMetrics.("UL Throughput Bytes")/simulationTime; ueName = [uesAssociatedTogNB.Name]'; gNBName = repmat(gnbs(m).Name,size(ueName,1),1); metricsSummarygNB = table(gNBName,ueName,dlTh,dlAvgSINR,dlMCS,dlPeakNumLayers,ulTh,ulAvgSINR,ulMCS,ulPeakNumLayers, ... VariableNames=["gNB" "UE" "DL Thput (Mbps)" "Avg DL SINR (dB)" "DL MCS" "Max DL NumLayers" "UL Thput (Mbps)" "Avg UL SINR (dB)" "UL MCS" "Max UL NumLayers"]); metricsSummary = [metricsSummary; metricsSummarygNB]; end end end function [connectionMatrix,gNBsTable,uesTable] = generateConnectionMatrix(gnbSites,ueSites,gNBsTable,uesTable) % Generate a connection matrix between gNBs and UEs. connectionMatrix is a % matrix of size numUEs by numgNBs, where numUEs is the number of UEs % (rxsites) in the ueSites array, and numgNBs is the number of gNBs % (txsites) in the gnbSites array. A 1 entry in the matrix indicates that % the corresponding gNB (row number) is connected to the UE corresponding % to the column number. % % A connection is established based on distance. UEs are connected to the % closest gNB. % Distance between all gNBs and all UEs d = distance(gnbSites,ueSites); [~,idxgnb] = min(d,[],2); connectionMatrix = ((1:size(d,2)) == idxgnb); % Add connection information to the gNB and UE tables for n=1:size(connectionMatrix,2) % for all gNBs gnbcol = connectionMatrix(:,n); if any(gnbcol) rowLogicalIdx = gNBsTable.gnbArrayIdx==n; gNBsTable.AssociatedUEsIdx(rowLogicalIdx) = {find(gnbcol)'}; gNBsTable.AssociatedUEs(rowLogicalIdx) = {uesTable.ueLabel(gnbcol)'}; end end % Connection matrix to table for n=1:size(connectionMatrix,1) % for all UEs uerow = connectionMatrix(n,:); if any(uerow) logicalIdx = uesTable.ueArrayIdx==n; uesTable.AssociatedgNBsIdx(logicalIdx) = {find(uerow)'}; uesTable.AssociatedgNBs(logicalIdx) = string(gnbSites(logical(uerow)).Name); end end end function [txUL,rxUL] = createULSites(tx,rx,PtxUEW,ulFc) % Turn tx sites into rxsites. This is used to perform ray tracing in the % uplink % Create rxsites for the gNBs ngNBs = numel(tx); rxUL = rxsite.empty(); for n=ngNBs:-1:1 rxArrayUL = tx(n).Antenna.clone(); rxUL(n) = rxsite(Latitude=tx(n).Latitude,Longitude=tx(n).Longitude, ... Antenna=rxArrayUL,AntennaAngle=tx(n).AntennaAngle, ... AntennaHeight=tx(n).AntennaHeight); end % Create txsites for the UEs nUEs = numel(rx); fc = ulFc; txUL = txsite.empty(); for n=nUEs:-1:1 txArrayUL = rx(n).Antenna.clone(); txUL(n) = txsite(Latitude=rx(n).Latitude,Longitude=rx(n).Longitude, ... TransmitterPower=PtxUEW,TransmitterFrequency=fc, ... Antenna=txArrayUL,AntennaHeight=rx(n).AntennaHeight,AntennaAngle=rx(n).AntennaAngle); end end function createRayTracingChannels(rtChannelRegistry,raysDL,raysUL,txgnbDL,rxueDL,txueUL,rxgNBUL,gnbs,ues) % Create a ray tracing channel for each gNB-UE pair. Skip channel creation % if no rays exist between the gNB and the UE. numgNBs = numel(gnbs); numUEs = numel(ues); % Create uplink and downlink ray tracing channels for gNBIdx = 1:numgNBs % Create channels between all UEs and all gNBs for ueIdx=1:numUEs % Downlink channels. Create a ray tracing channel model per UE if ~isempty(raysDL{gNBIdx}{ueIdx}) rtChannelRegistry.addRayTracingChannel(raysDL{gNBIdx}{ueIdx},txgnbDL(gNBIdx),rxueDL(ueIdx),gnbs(gNBIdx).ID,ues(ueIdx).ID); end % Uplink channels if ~isempty(raysUL{gNBIdx}{ueIdx}) rtChannelRegistry.addRayTracingChannel(raysUL{gNBIdx}{ueIdx},txueUL(ueIdx),rxgNBUL(gNBIdx),ues(ueIdx).ID,gnbs(gNBIdx).ID); end end end end function [mcs,sinr,numLayers,slotNum] = getMCSSINR(ueLogs,slotsPerFrame) % Extract the MCS, SINR and absolute slot number from the UE logs UELOGS % given the number of slots per frame SLOTSPERFRAME. % Calculate absolute slot number fn = ueLogs.Frame; sn = ueLogs.Slot; slotN = fn.*slotsPerFrame + sn; % Sort absolute slot number in ascending order [slotNum,sortIdx] = sort(slotN); % Extract and reorder MCS and SINR mcs = ueLogs.("MCS Index"); sinr = ueLogs.SINR; numLayers = ueLogs.NumLayers; mcs = mcs(sortIdx); sinr = sinr(sortIdx); end function raysUL = dl2ulRays(rays) % Convert downlink rays to uplink rays when we assume channel reciprocity. % This function swaps the transmitter and receiver locations, "flips" the % interaction matrix, and transposes the polarization matrix. Other uplink % ray properties remain the same. Constructing the uplink rays this way % results on the swapping of the angles of arrival and departure. raysUL = rays; % Allocate for n = 1:numel(rays) % loop over all receivers for this transmitter for m = 1:numel(rays{n}) % loop over all rays for this tx/rx link raysUL{n}(m).TransmitterLocation = rays{n}(m).ReceiverLocation; raysUL{n}(m).ReceiverLocation = rays{n}(m).TransmitterLocation; raysUL{n}(m).PolarizationMatrix = rays{n}(m).PolarizationMatrix.'; raysUL{n}(m).Interactions = flip(rays{n}(m).Interactions); end end end
