-
Notifications
You must be signed in to change notification settings - Fork 12
/
Copy pathTheoryLOS.m
304 lines (237 loc) · 11.4 KB
/
TheoryLOS.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
% Written by Ish Jain
% NYU Tandon School of Engineering
% Date: June 2018
%
% Description:
%
% Use generalized model and find dynamic blockage probability given
% coverage. Coverage is defined as "atleast one BS out side blocked zone
% (defined by static and self-blockage)"
clear
close all
wannaplotCellRadius=0; %change only when focussing on cell radius
wannaSaveFiles=1; %to save csv files for plotting using Latex.
wannaplot=1; %plot results to visualize
V = 1; %velocity of blocker m/s
hb = 1.8; %height of blocker
hr = 1.4; % height UE
ht = 5;%height BS
frac = (hb-hr)/(ht-hr); %fraction depends on heights
mu = 2; %Expected bloc dur =1/mu
densityBS = [50,100,200,300,400,500]*10^(-6); %BS=Base Station
densityBL = [0.01, 0.1]; %Dynamic BLockers
densityD = [1e-9,0.0001]; %D = static blockage
omegaVal = [0 pi/3]; %self-blockage angle
nBS = length(densityBS);
nBL = length(densityBL);
nD = length(densityD);
nO = length(omegaVal);
El = 10; %building length %m
Ew = 10; %building width %m
pB = zeros(nBS,nBL,nO);
pBCond=zeros(nBS,nBL,nO);
freqCond=zeros(nBS,nBL,nO);
durCond=zeros(nBS,nBL,nO);
R = 200; %m Radius
for iT = 1:nBS
tempind = 0;
for iB = 1:nBL
for iD = 1:nD
for iO = 1:nO
tempind=tempind+1; %increment from zero after every BSdensity value
lamT = densityBS(iT); %lambda BS
lamB = densityBL(iB);
lamD = densityD(iD);
omega = omegaVal(iO);
p(iO) =1- omega/(2*pi);
C(iB) = 2/pi*lamB*V*frac;
a(iB) = 2*mu/(R*C(iB))-2*mu^2/(R^2*C(iB)^2)*log(1+R*C(iB)/mu);
% Coverage probability for open park like area
pnn0(iT,iO)=exp(-p(iO)*lamT*pi*R^2); %prob that n not equal to 0 (pnn0)
% Define parameters for static blockage
beta(iD) = 2/pi*lamD*(El+Ew);
beta0(iD) = lamD*El*Ew;
b(iD) = 2*exp(-beta0(iD))/(beta(iD)^2*R^2) * ...
(1-(1+beta(iD)*R)*exp(-beta(iD)*R));
%write code to verify b: verified. Both matches
bExp = @(r) exp(-(beta(iD).*r+beta0(iD))) *2.*r/R^2;
b(iD) = integral(bExp,0,R);
%coverage probability given by static and self blockage
pCoverage(iT,iD,iO) = 1-exp(-b(iD)*p(iO)*lamT*pi*R^2);
pNoCoverage(iT,iD,iO) = 1-pCoverage(iT,iD,iO);
%1. Marginal prob of Blockage (open area)
pB(iT,iB,iO) = exp(-a(iB)*p(iO)*lamT*pi*R^2);
%1.1 considering static blockage too.
c=C(iB)/mu;
aPrimeExp = @(r) exp(-(beta(iD)*r+beta0(iD)))./(1+c*r) *2.*r/R^2;
aPrime(iB,iD) = integral(aPrimeExp, 0, R);
pBprime(iT,tempind) = exp(-aPrime(iB,iD)*p(iO)*lamT*pi*R^2);
%2. Conditional prob of blockage given coverage(newly defined)
pBCondPrime(iT,tempind) = (pBprime(iT,tempind)-...
pNoCoverage(iT,iD,iO))/pCoverage(iT,iD,iO);
%3. Expected Freq of blockage
freqPrime(iT,tempind) = mu*(1-aPrime(iB,iD))*p(iO)*b(iD)...
*lamT*pi*R^2*exp(-aPrime(iB,iD)*p(iO)*b(iD)*lamT*pi*R^2);
%4. Conditional expectation of freq of bl given n!=0
freqCondPrime(iT,tempind) = freqPrime(iT,tempind)/pCoverage(iT,iD,iO);
%5. Conditional expectation of duration of bl given n!=0
%This is used to get duration
tempb(iT,iD,iO) = p(iO)*b(iD)*lamT*pi*R^2;%*angFrac;
Ei(iT,iB,iO) = ei(tempb(iT,iD,iO))-log(tempb(iT,iD,iO))-0.5772; %exponential integral function
durCondPrime(iT,tempind) = (1-pCoverage(iT,iD,iO))*Ei(iT,iB,iO)/(mu*pCoverage(iT,iD,iO));
durCondPrime(iT,tempind) = durCondPrime(iT,tempind)*1000;
%Try different approximations
durCond(iT,iB,iO) = (1/tempb(iT,iD,iO)+1/tempb(iT,iD,iO)^2)/...
(mu*(1-pnn0(iT,iO)));
durCond2(iT,iB,iO) = (1/tempb(iT,iD,iO))/...
(mu*(1-pnn0(iT,iO)));
durCond2Prime(iT,tempind) = (1/tempb(iT,iD,iO))/...
(mu*pCoverage(iT,iD,iO));
durCond2Prime(iT,tempind)=durCond2Prime(iT,tempind)*1000;
colTitle{1}='lamT';
if(lamD<1e-7),lamD=0;end
colTitle{tempind+1} = strcat('lamB',num2str(lamB),...
'lamD',num2str(lamD*1e4),'omega',num2str(omega*360/2/pi));
legendArray{tempind} = strcat(' \lambda_B=',num2str(lamB),...
'\lambda_D=',num2str(lamD*1e4),' \omega=',num2str(omega*360/2/pi));
% sprintf('\lambda_B=%f,\lambda_D=%d,\omega=%d',lamB,lamD*1e6,omega*360/2/pi);
end
end
end
end
if(wannaSaveFiles)
writetable(cell2table([colTitle; num2cell([densityBS'*10^4,pBprime])]),...
'figures/theory_pB.csv','writevariablenames',0)
writetable(cell2table([colTitle; num2cell([densityBS'*10^4,pBCondPrime])]),...
'figures/theory_pBCond.csv','writevariablenames',0)
writetable(cell2table([colTitle; num2cell([densityBS'*10^4,freqPrime])]),...
'figures/theory_freq.csv','writevariablenames',0)
writetable(cell2table([colTitle; num2cell([densityBS'*10^4,freqCondPrime])]),...
'figures/theory_freqCond.csv','writevariablenames',0)
writetable(cell2table([colTitle; num2cell([densityBS'*10^4,durCond2Prime])]),...
'figures/theory_durCond.csv','writevariablenames',0)
end
if(wannaplot)
figure(1);grid on;
semilogy(densityBS,pBprime);
ylim([1e-6,1]);title('Marginal prob of Blockage')
legend(legendArray);
figure(2);grid on;
semilogy(densityBS,pBCondPrime);
title('Conditional prob of Bl given n!=0');
ylim([1e-6,1])
legend(legendArray);
figure(3);grid on;
semilogy(densityBS,freqPrime)
title('Expected Freq of blockage')
legend(legendArray);
ylim([1e-6,1])
figure(4);grid on;
semilogy(densityBS,freqCondPrime);
title('Conditional expectation of freq of bl given n!=0')
ylim([1e-6,1])
legend(legendArray);
figure(5); grid on;
plot(densityBS,durCond2Prime)
hold on
% plot(densityAP',durCond2(:,6),'y-')
% plot(densityAP',durCond3(:,6),'b-')
title('Conditional expectation of duration of bl given n!=0')
legend(legendArray);
end
%% Effect of cell radius
%Just copied the above code with minor modifications
if(wannaplotCellRadius)
clear
close all
wannaplot=1;
V = 1; %velocity m/s
hb = 1.8;
hr = 1.4;
ht = 5;
frac = (hb-hr)/(ht-hr); %fraction depends on heights
mu = 2; %Expected bloc dur =1/mu
R = 100; %m Radius
densityBS = [0.01,1,50,100,150,200,250,300,350,400]*10^(-6); %BS=Base Station
% densityBS = [100]*10^(-6); %BS=Base Station
densityBL = [ 0.01]; %Dynamic BLockers
densityD = [0.0001]; %D = static blockage
omegaVal = [ pi/3];
% Rvalues = 50:500; %m Radius
Rvalues = [100, 200];
nR = length(Rvalues);
nBS = length(densityBS);
nBL = length(densityBL);
nD = length(densityD);
nO = length(omegaVal);
El = 10; %m
Ew = 10; %m
% !!!-------------------------------!!!
pB = zeros(nBS,nBL,nO);
pBCond=zeros(nBS,nBL,nO);
freqCond=zeros(nBS,nBL,nO);
durCond=zeros(nBS,nBL,nO);
for iT = 1:nBS
tempind = 0;
for iR = 1:nR
R = Rvalues(iR);
for iB = 1:nBL
for iD = 1:nD
for iO = 1:nO
tempind=tempind+1; %increment from zero after every BSdensity value
lamT = densityBS(iT); %lambda BS
lamB = densityBL(iB);
lamD = densityD(iD);
omega = omegaVal(iO);
p(iO) =1- omega/(2*pi);
C(iB) = 2/pi*lamB*V*frac;
a(iB) = 2*mu/(R*C(iB))-2*mu^2/(R^2*C(iB)^2)*log(1+R*C(iB)/mu);
% Coverage probability for open park like area
pnn0(iT,iO)=exp(-p(iO)*lamT*pi*R^2);
beta(iD) = 2/pi*lamD*(El+Ew);
beta0(iD) = lamD*El*Ew;
b(iD) = 2*exp(-beta0(iD))/(beta(iD)^2*R^2) * ...
(1-(1+beta(iD)*R)*exp(-beta(iD)*R));
%write code to verify b: verified. Both matches
bExp = @(r) exp(-(beta(iD).*r+beta0(iD))) *2.*r/R^2;
b(iD) = integral(bExp,0,R);
%coverage probability given by static and self blockage
pCoverage(iT,iD,iO) = 1-exp(-b(iD)*p(iO)*lamT*pi*R^2);
pNoCoverage(iT,iD,iO) = 1-pCoverage(iT,iD,iO);
%1. Marginal prob of Blockage (open area)
pB(iT,iB,iO) = exp(-a(iB)*p(iO)*lamT*pi*R^2);
%1.1 considering static blockage too.
c=C(iB)/mu;
aPrimeExp = @(r) exp(-(beta(iD)*r+beta0(iD)))./(1+c*r) *2.*r/R^2;
aPrime(iB,iD) = integral(aPrimeExp, 0, R);
pBprime(iT,tempind) = exp(-aPrime(iB,iD)*p(iO)*lamT*pi*R^2);
%2. Conditional prob of Bl given n!=0
pBCond(iT,iB,iO) = (pB(iT,iB,iO)-pnn0(iT,iO))./(1-pnn0(iT,iO));
%2.1 Conditional prob of blockage given coverage(newly defined)
pBCondPrime(iT,tempind) = (pBprime(iT,tempind)-...
pNoCoverage(iT,iD,iO))/pCoverage(iT,iD,iO);
colTitle{1}='Radius';
if(lamD<1e-7),lamD=0;end
colTitle{tempind+1} = strcat('lamB',num2str(lamB),...
'lamD',num2str(lamD*1e4),'omega',num2str(omega*360/2/pi),...
'R',num2str(R));
%Put legends
legendArray{tempind} = strcat(' \lambda_B=',num2str(lamB),...
'\lambda_D=',num2str(lamD*1e4),' \omega=',num2str(omega*360/2/pi),...
'R=',num2str(R));
end
end
end
end
end
% writetable(cell2table([colTitle; num2cell([Rvalues', pBCondPrime])]),...
% 'figures/theory_withR_LOS.csv','writevariablenames',0);
% save('figures/R_LOS.mat','pBCondPrime');
figure(8);grid on;
semilogy(densityBS,pBCondPrime, 'LineWidth',2);
xlabel('BS Density');
title('LOS Blockage Probability');
ylim([1e-6,1])
legend(legendArray);
% legend('R=100m','R=200m','R=500m','R=1000m');
end