Simulating Cochannel Interference in Wireless Networks

I did a project recently that involved the simulation of cochannel interference in a cellular wireless network that I found pretty interesting. It was a monte carlo simulation that randomly positioned a mobile device within a hexagonal cell (shown in green) inside a network with six other cochannels.

The goal was to calculate the signal-to-interference (SI) ratio of the mobile device and calculate the probability of a user experiencing degraded service. The measure of degraded service was defined as the probability of the SI ratio being below 18dB. The propogation environment environment was characterized by a log-distance path loss model with shadowing.

After running the simulation for several thousand iterations of varying path loss exponents and noise levels I obtained plots of the probability distribution of the SI ratio. The green vertical line is the threshold of 18dB and the blue line is the mean of the SI ratio.

Here's the code used to generate the results.

Matlab Code

close all
clear all
clcformat short g

num_samples = 16384; % more is better
N = 7; % Freq. reuse factor
Ni = 2; Nj = 1; % N = i^2+i*j+j^2
R = 10000; % Cell radius in m
SI_thresh = 18; % in dB
Pt = 100;  % Transmitter power in dBW
Lpd0 = 5; % Path loss at reference d0 in dBW
d0 = 500; % ref distance in m
d00 = 1/d0; % divides are time consuming so use this in the simulation

se = [7 8 9];   % sigma_epsilon in dB
num_se = length(se);
k = [2, 2.5, 3, 3.5, 4];
num_k = length(k);
num_cells = 6;
% Calculate the location of cochannel cells
theta = pi/180*[30 90 150 -150 -90 -30].';
result = sqrt(3)*R*(Ni*exp(1i*theta) + Nj*exp(1i*(theta +60*pi/180)));
cell_pos = [real(result) imag(result)];
% Plot cells
if 1
 % Display a plot of the cell with inner circle
 hex_theta = 0:60:360;
 plot(R*cosd(hex_theta), R*sind(hex_theta))
 hold on
 cp = d0*exp(1i*linspace(0, 2*pi));
 plot(real(cp), imag(cp))
 axis('square')

 % Plot the cochannel hexagons
 if 1
   for idx=1:size(cell_pos,1)
     x = cell_pos(idx, 1);
     y = cell_pos(idx, 2);
     plot(R*cosd(hex_theta)+x, R*sind(hex_theta)+y)
     plot(x, y, 'xk')
   end
 end
end
% Preallocate some matrices (for speed) (I should preload a randn matrix,
% but that would add too much complexity to the code)
samples = zeros(length(k), num_samples);
positions = zeros(num_samples, 2);
bad_positions = zeros(num_samples, 2);
dcells = zeros(num_cells, 1);

profile on
tic% Start the simulation
idx = 1;
idx2 = 1;
while(idx <= num_samples)
 % Generate the position of a MH within the hexagonal cell
 if 1
   % part a; position is random inside the cell
   xpos = -R*(1 - 2*rand(1)); % [+/-R]
   ypos = -0.5*sqrt(3)*R*(1 - 2*rand(1)); % [+/-sqrt(3)/2*R]
 else
   % Part b: position is at the edge of the cell
   xpos = R;
   ypos = 0;
 end
   
 % Determine if this is a valid position
 % d > do and inside the area of the hexagon
 d = sqrt(xpos^2 + ypos^2);
 y = -sqrt(3)*abs(xpos) + sqrt(3)*R;

 if (d < d0 || abs(ypos) > y)
   % Save this for plotting purposes
   bad_positions(idx2, 1) = xpos;
   bad_positions(idx2, 2) = ypos;
   idx2 = idx2 + 1;
   continue;
 else
   % Save this for plotting purposes
   positions(idx, 1) = xpos;
   positions(idx, 2) = ypos;
 end

 % Calculate the distance between cochannels and the MH
 dcells = sqrt((cell_pos(:,1)-xpos).^2 + (cell_pos(:,2)-ypos).^2);

 % Loop over all values of path loss (k) and noise level (se) calculating
 % the signal-to-interference level in dB
 for kidx = 1:num_k
   for seidx = 1:num_se
     % Calculate S
     % Lp(d) = L_bar_p(d0) + 10*k*log10(d/d0) + noise (dB)
     S = Pt - (Lpd0 + 10*k(kidx)*log10(d*d00) + se(seidx)*randn(1));

     % Calculate I (sum of I's))
     % Lp(d) = L_bar_p(d0) + 10*k*log10(d/d0) + noise (dB)
     I = Pt - (Lpd0 + 10*k(kidx)*log10(dcells*d00) + se(seidx)*randn(6,1));
     I = 10*log10(sum(10.^(0.1*I)));
 
     % S/I = Signal - Interference (in dB)
     samples(kidx, idx, seidx) = S - I;
   end
 end

 % Increment the storage/simulation index
 idx = idx + 1;
end
toc
profview
profile off
% Plot the good and bad positions for MH
plot(positions(:, 1), positions(:, 2), 'gx')
plot(bad_positions(:, 1), bad_positions(:, 2), 'rx')
% Print out some statistics of the simulations as well as P[S/I < Thresh]
result = squeeze(sum(samples < SI_thresh, 2)) ./ num_samples;
result = [k(:) result];
result = [[0 se]; result];
disp(result)
% Subplot of all the distributions resulting from the simulation
xmin = min(samples(:));
xmax = max(samples(:));
ymin = 0; ymax = 0;
figure
for kidx=1:num_k
 for sidx=1:num_se
   % Get the current subplot
   subplot(num_k, num_se, (kidx-1)*num_se + sidx)
   
   % Compute the histogram based pdf
   [n, x] = hist(samples(kidx, :, sidx), 32);
   n = n./trapz(x, n);
   
   % Compute the mean so I can display a vertical reference line for plots
   mean_sample = mean(samples(kidx, :, sidx));
   
   % plot the pdf, threshold, and the mean line
   plot(x, n)
   line([mean_sample mean_sample], [min(n) max(n)], 'color', 'b', 'LineWidth', 2)
   line([SI_thresh SI_thresh], [min(n) max(n)], 'color', 'g', 'LineWidth', 2)
       
   title(sprintf('k: %.1f sigma: %.1f', k(kidx), se(sidx)))
   grid on, hold on
   
   % Compute ymin/ymax so I can make all the axis limits equal (later)
   ymin = min(min(n), ymin); ymax = max(max(n), ymax);
 end
end
% This code resets the subplot axis from the previous loop
for kidx=1:num_k
 for sidx=1:num_se
     subplot(num_k, num_se, (kidx-1)*num_se + sidx)
     axis([xmin xmax ymin ymax*1.1])
 end
end