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
| function main
close all;
%% Configurate parameters
nTrials = 500;% Number of Trials
nT = 2; % N Tx antennas
nR = 2; % N Rx antennas
nEveArray = [2,3,4,5,6,7]; % N Eve antennas
% SNR array
snrArray_db = [-10, 10]; % dB
snrArray_lin = db2lin(snrArray_db); % Power in mW
nSnr = numel(snrArray_db); % N SNR points
% Seed
rng(1);
% Step for lymbda optimization
stepForLymbda = 0.005;
%% Generate lymbda array
lymbdaArray = 0:stepForLymbda:1;
nLyambdaMax = numel(lymbdaArray);
lyambda1Arr = ((0:stepForLymbda:1)' * ones(1,nLyambdaMax))';
lyambda2Arr = ((1:-stepForLymbda:0)' * ones(1,nLyambdaMax))';
%% Allocate arrays
rArray = zeros(numel(snrArray_db), nTrials);
lymbda1 = zeros(nLyambdaMax,nLyambdaMax,nSnr);
lymbda2 = zeros(nLyambdaMax,nLyambdaMax,nSnr);
%% Generate lymbda arrays for each SNR
for iSnr = 1 : numel(snrArray_db)
snr_lin = snrArray_lin(iSnr);
lymbda1(:,:,iSnr) = lyambda1Arr * snr_lin;
lymbda2(:,:,iSnr) = lyambda2Arr * snr_lin;
end
%% Main loop
for iE = 1 : numel(nEveArray)
nE = nEveArray(iE)
for iTrial = 1:nTrials
%iTrial
% Generate random channels
H = genChannel(nT,nR);
G = genChannel(nT,nE);
% Find channel matrix values
HhT = H*H';
GgT = G*G';
h1 = HhT(1,1);
h2 = HhT(1,2);
h3 = HhT(2,2);
g1 = GgT(1,1);
g2 = GgT(1,2);
g3 = GgT(2,2);
% Find additional parameters
a1 = (lymbda2 - lymbda1)./h2;
b1 = 1/2.*(lymbda1 - lymbda2)./(h3-h1);
c1 = 1 + 1/2*(lymbda1 + lymbda2).*(h1+h3)+lymbda1.*lymbda2.*(h1*h3 - h2*h2);
a2 = (lymbda2 - lymbda1).*g2;
b2 = 1/2.*(lymbda1 - lymbda2).*(g3-g1);
c2 = 1 + 1/2.*(lymbda1+lymbda2).*(g1+g3) + lymbda1.*lymbda2.*(g1*g3-g2*g2);
a = c1.*b2 - c2.*b1;
b = a1.*c2 - a2.*c1;
c = a1.*b2 - a2.*b1;
% Find optimal tetta
tetta = -1/2*atan(b./a) + 1/2*asin(c./sqrt(a.*a+b.*b)) + pi/2;
% Calculate r
realValR = ((((a1.*sin(2.*tetta)+b1.*cos(2.*tetta) + c1)./(a2.*sin(tetta.*2)+b2.*cos(tetta.*2) + c2))));
% Maximize r per lymbdas
maxRealValR = max(realValR,[],1);
maxRealValR = squeeze( max(maxRealValR,[],2));
% Save r to the array
rArray(:,iTrial) = maxRealValR;
end
% Calculate r and find mean value
r = mean(1/2.*log2(rArray),2);
%% Plot figures
% Plot r per SNR
figure;
plot(snrArray_db,r,'-*', 'LineWidth',2);
ylim([min(0,min(r)),max(r)]);
grid on;
legendline = (strcat('nTx = ', num2str(nT), ', nRx = ', num2str(nR), ', nE =', num2str(nE)));
legend(legendline);
xlabel('SNR, dB')
ylabel('R, bps/Hz')
title('R per SNR')
savefig(strcat('R per SNR','_',legendline));
% Plot r per Power
% figure;
% plot(snrArray_lin,r,'-*', 'LineWidth',2);
% ylim([min(0,min(r)),max(r)]);
% grid on;
% legendline = (strcat('nTx = ', num2str(nT), ', nRx = ', num2str(nR), ', nE =', num2str(nE)));
% legend(legendline);
% xlabel('Power, mW')
% ylabel('R, bps/Hz')
% title('R per Tx power in case n_e = 1 mW')
% savefig(strcat('R per Tx power in case ne = 1 mW','_',legendline));
end
end
function val_lin = db2lin(val_db)
val_lin = 10.^(val_db./10);
end
function channelMatrix = genChannel(nT,nR)
channelMatrix = randn(nT,nR);
end |