Analyze SINR in 6G FR3 Network Using Digital Twin

nUEs = 20;              % Number of UEs
fc = 14e9;              % Carrier frequency
PtxdBm = 40;            % gNB transmit power in dBm
PtxInterfdBm = 43;      % Interfering gNB transmit power in dBm
Noc = -145;             % dBm/Hz from 38.101
NRB = 52;               % Number of resource blocks
maxNumReflections = 3;  % Maximum number of ray tracing reflections
maxNumDiffractions = 1; % Maximum number of ray tracing diffractions
txArraySize = [4 4];    % Number of antennas (rows, columns) in transmitter (gNB) rectangular array
rxArraySize = [2 2];    % Number of antennas (rows, columns) in receiver (UE) rectangular array

rng("default");         % Random number generator default settings for result reproducibility

filename = "chicago.osm";

latROILim = [41.8845 41.8781];
lonROILim = [-87.6354 -87.6279];

ROI = RegionOfInterest(filename,latROILim,lonROILim);

ueLocations = ROI.GenerateRandomDrops(nUEs);

% Plot generated UE location points on map
ROI.plotMap(true,ueLocations)

if exist("viewer","var") && isvalid(viewer) % Viewer handle exists and viewer window is open
    clearMap(viewer);
else
    viewer = siteviewer(Basemap="openstreetmap",Buildings=filename);
end

gNBLatLon = [41.881696 -87.629566];
lambda = physconst("LightSpeed")/fc;
txArray = phased.URA(Size=txArraySize,ElementSpacing=0.5*lambda);
PtxW = db2pow(PtxdBm-30);
tx = txsite(Latitude=gNBLatLon(1),Longitude=gNBLatLon(2), ...
    TransmitterPower=PtxW, ...
    TransmitterFrequency=fc, ...
    Antenna=txArray, ...
    AntennaAngle=90);

show(tx)

% Create all rxsites
rx = rxsite(Latitude=ueLocations(:,1),Longitude=ueLocations(:,2));
% Assign arrays to each rxsite
for n = 1:nUEs
    rx(n).Antenna = phased.URA(Size=rxArraySize,ElementSpacing=0.5*lambda);
end

pm = propagationModel("raytracing",MaxNumReflections=maxNumReflections,MaxNumDiffractions=maxNumDiffractions);

rays = raytrace(tx,rx,pm);

% Carrier object to obtain sampling rate and OFDM FFT size
carrier = pre6GCarrierConfig(NSizeGrid=NRB);

% Calculate channel matrices using SVD
H = calculateChannelMatrix(carrier,tx,rx,rays);

% Calculate precoding (digital beamforming) and combining weights for all gNB to UE links
[Wrx,S,Wtx] = pagesvd(H);

% UE received signal power (dBm)
ueSignalPowerdBm = calculateSignalPower(tx,rx,rays,Wtx,Wrx);

% UE received signal power (W)
ueSignalPowerW = dBm2W(ueSignalPowerdBm);

interfBS = [41.884085,-87.637310;
            41.885208,-87.625557;
            41.877766,-87.635363;
            41.877684,-87.625858];
txInterf = txsite(Latitude=interfBS(:,1),Longitude=interfBS(:,2), ...
    TransmitterFrequency=fc, ...
    TransmitterPower=dBm2W(PtxInterfdBm));

show(txInterf)

% Assign MIMO combining weights to each UE before calculating the received
% interfering power
isLoS = zeros(nUEs,1);
for ueIdx=1:nUEs
    wrx = Wrx(:,1,ueIdx).'; % Weights from strongest singular value
    rx(ueIdx).Antenna.Taper = conj(wrx);
    isLoS(ueIdx) = any([rays{ueIdx}.LineOfSight]);
end

% Calculate interfering power for all UEs from each interfering gNB
numInterf = length(txInterf);
interfPowerdBm = zeros(nUEs,numInterf);
for n=1:numInterf
    interfPowerdBm(:,n) = sigstrength(rx,txInterf(n),pm);
end

% Interference power per UE per interfering base station
interfPowerW = dBm2W(interfPowerdBm);

% Combine the interference power for all interfering base stations per UE
% (in W)
ueInterfPowerW = sum(interfPowerW,2);

% Noise power
noisePowerdBm = Noc+10*log10(carrier.NSizeGrid*12*carrier.SubcarrierSpacing*1e3);
noisePowerW = dBm2W(noisePowerdBm);

% Calculate SINR
SINR = 10*log10(ueSignalPowerW./(ueInterfPowerW+noisePowerW));

% Result table
varNames = ["UE" "UE signal power (dBm)" "LOS" "Interf power (dBm)" "Noise+Interf (dBm)" "SINR"];
SINRTable = table((1:nUEs).',ueSignalPowerdBm,isLoS,W2dBm(ueInterfPowerW),W2dBm(ueInterfPowerW+noisePowerW),SINR,VariableNames=varNames);
disp(SINRTable);

ROI.plotMap(true,ueLocations,SINR,[tx.Latitude tx.Longitude],interfBS)

function yW = dBm2W(ydBm)
% dBm to W conversion
    yW = 10.^((ydBm-30)/10);
end

function ydBm = W2dBm(yW)
% W to dBm conversion
    ydBm = pow2db(yW)+30;
end

function H = calculateChannelMatrix(carrier,tx,rx,rays)
% Calculate the channel matrices for all links between the gNB and all
% receivers using the result of the ray tracing analysis. The resulting
% channel matrix H is of size nRx-by-nTx-by-nUEs, where nRx and nTx are the
% number of transmit and receive antennas respectively, and nUEs is the
% number of UEs. The channel matrix is obtained using an average in time
% and frequency over a slot.

    % Obtain sampling rate and OFDM FFT size
    ofdmInfo = nrOFDMInfo(carrier);
    fs = ofdmInfo.SampleRate;
    
    % Number of transmit and receive antennas
    nTx = prod(tx.Antenna.Size);
    nRx = prod(rx(1).Antenna.Size); % Assume all receivers have the same number of antennas
    
    % Tx power per antenna
    PtxPerAnt = tx.TransmitterPower/nTx;
    
    % Calculate channel matrix
    nUEs = length(rx);
    H = zeros(nRx,nTx,nUEs);
    for ueIdx=1:nUEs
        if ~isempty(rays{ueIdx})
            % Create channel object
            ch = comm.RayTracingChannel(rays{ueIdx},tx,rx(ueIdx));
            ch.SampleRate = fs;
            ch.NormalizeImpulseResponses = false; % maintain path loss information
            ch.NormalizeChannelOutputs = false;
            ch.ChannelFiltering = false;
            ch.NumSamples = sum(ofdmInfo.SymbolLengths(1:14)); % one slot
            cir = ch(); % channel impulse response
            chInfo = ch.info();
    
            % OFDM channel response
            pathFilters = chInfo.ChannelFilterCoefficients;     
            K = ofdmInfo.Nfft.^2/(carrier.NSizeGrid*12);
            hest = sqrt(nTx)*sqrt(PtxPerAnt)*sqrt(K)*nrPerfectChannelEstimate(carrier,cir,pathFilters.'); 
    
            % Average channel response over time and freq
            H(:,:,ueIdx) = squeeze(mean(hest(:,:,:,:),[1 2]));
        end
    end
end

function ueSignalPowerdBm = calculateSignalPower(tx,rx,rays,Wtx,Wrx)
% Calculate the received signal power given a txsite, nUEs rxsites
% (rx), using the rays between the transmitter and receivers. Wtx and Wrx
% are the precoding and combining matrices to consider for the power
% calculation.

    nUEs = length(rx);
    fc = tx.TransmitterFrequency;
    PtxdBm = W2dBm(tx.TransmitterPower);
    ueSignalPowerdBm = zeros(nUEs,1);
    for ueIdx=1:nUEs
        if ~isempty(rays{ueIdx})
            % Precoding and combining weights corresponding to this gNB to
            % UE link
            wtx = Wtx(:,1,ueIdx).'; % Weights from strongest singular value
            wrx = Wrx(:,1,ueIdx).'; % Weights from strongest singular value
    
            % AoD and AoA angles
            AoD = [rays{ueIdx}.AngleOfDeparture];
            AoD(1,:) = AoD(1,:) - tx.AntennaAngle;
            AoA = [rays{ueIdx}.AngleOfArrival];
            AoA(1,:) = AoA(1,:) - rx(ueIdx).AntennaAngle;
    
            % Make sure rx array has no taper, as weights are specified as
            % an argument when calling the directivity() function, and not
            % as part of the array taper property.
            rx(ueIdx).Antenna.Taper = 1;
    
            % Calculate the received signal power per UE
            txDirectivity = directivity(tx.Antenna, fc, wrapTo180(AoD),"Weights",wtx.');
            rxDirectivity = directivity(rx(ueIdx).Antenna, fc, wrapTo180(AoA),"Weights",wrx.');
            rxPowerdBW = PtxdBm - 30 + txDirectivity' + rxDirectivity' - [rays{ueIdx}.PathLoss]; % in dBW
            rxSignalLevel = exp(1i*[rays{ueIdx}.PhaseShift]) .* 10.^(rxPowerdBW/20); % in V
            ueSignalPowerdBm(ueIdx) = W2dBm(abs(sum(rxSignalLevel)).^2); % dBm
        else
            ueSignalPowerdBm(ueIdx) = -Inf;
        end
    end
end