仿真内容
文件中包含了一个关于无源多用户接入（Unsourced Multiple Access，UMA）系统的 MATLAB 数值例程，用于评估随机用户活动情况下的随机编码界限。
这个工作主要在论文 [1] 中介绍，该论文题为 "Unsourced Multiple Access With Random User Activity"，于 2022 年 1 月提交给 IEEE Transactions on Information Theory。
这个工作的部分内容也在论文 [2] 中介绍，该论文题为 "Massive Uncoordinated Access With Random User Activity"，于 2021 年 7 月在 IEEE International Symposium on Information Theory (ISIT) 中展示。
如果你使用了这些代码，请引用上述的论文。
代码介绍部分
/RCU_KaKnown/这个文件夹包含了关于无源多用户接入系统中的每用户误差概率（PUPE）的随机编码界限。具体包括以下文件：
RCU_KaFixedKnown.m对于活跃用户数 Ka 固定且已知的情况下，计算 PUPE 的随机编码界限。
RCU_KaRandomKnown.m对于活跃用户数 Ka 是随机的情况下，计算 PUPE 的随机编码界限。
RCU_KaPoissonKnown.m对于活跃用户数 Ka 服从泊松分布且已知的情况下，计算 PUPE 的随机编码界限。
EbN0_KaPoissonKnown.m找到使得活跃用户数 Ka 服从泊松分布且已知时 PUPE 的随机编码界限低于某个阈值的最小所需 EbN0。
binary_search.m,
golden_search.m辅助函数。
/RCU_KaUnknown/: 这个文件夹包含了关于无源多用户接入系统中，活跃用户数 Ka 是随机且未知的情况下的随机编码界限，这是论文 [1] 和 [2] 中提出的。具体包括以下文件：
RCU_KaRandomUnknown.m对于活跃用户数 Ka 是随机的情况下，计算 PUPE 的随机编码界限。
RCU_KaPoissonUnknown.m对于活跃用户数 Ka 服从泊松分布的情况下，计算 PUPE 的随机编码界限。
RCU_floor_KaRandomUnknown.m论文 [1, Theorem 3] 中对误差底线的特征化。
EbN0_KaPoissonUnknown.m找到使得活跃用户数 Ka 服从泊松分布且未知时 PUPE 的随机编码界限在某个阈值以下的最小所需 EbN0。
binary_search_P_MDFA.m,
golden_search_P1_MDFA辅助函数。
/SA-MPR/这个文件夹包含了上述界限在带有多包接收（MPR）的时隙ALOHA（SA-MPR）系统中的应用。
EbN0_SAMPR_KaPoissonKnown.m找到使得 SA-MPR 中活跃用户数 Ka 服从泊松分布且已知时的最小所需 EbN0。
EbN0_SAMPR_KaPoissonUnknown.m找到使得 SA-MPR 中活跃用户数 Ka 服从泊松分布且未知时的最小所需 EbN0。参考论文 [1, Corollary 2]。
/TIN/这个文件夹包含了一个将干扰视为噪声（Treat Interference as Noise，TIN）的方案的随机编码界限。
EbN0_TIN_KaPoissonUnknown.m找到使得 TIN 中活跃用户数 Ka 服从泊松分布且未知时的最小所需 EbN0。
注意: SPARC 和 enhance SPARC 方案的代码由 Krishna R. Narayanan 的研究组提供，因此未包含在此存储库中。
RCU_KaKnown
binary_search.m
function [epsilon, P] = binary_search(f,x1,x2,TARGET,TOL)
% search the value of P such that min_{P1} f(P,P1) is from TARGET - TOL to TARGET
iter = 20; % maximum number of iterations
k1 = 0; % iteratin index
% we use a golden search to compute min_{P1} f(P,P1)
fx = golden_search(f,1e-9,(x1+x2)/2,(x1+x2)/200);
while (TARGET < fx || fx<TARGET - TOL) && (k1<iter || TARGET < fx)
if k1 > 0
fx = golden_search(f,0,(x1+x2)/2,(x1+x2)/200);
end
if TARGET > fx
x2 = (x1+x2)/2; %set new end of interval
else
x1 = (x1+x2)/2; %replace as new start index
end
k1=k1+1;
end
epsilon = fx;
P = x1;
end
golden_search.m
function [epsilon, P1] = golden_search(f, START_INT, END_INT, TOL)
% minimize f(P1) over P1. epsilon = min_{P1} f(P1)
P = END_INT;
iter = 20; % maximum number of iterations
tau = double((sqrt(5)-1)/2); % golden proportion coefficient, around 0.618
k2 = 0; % iteration index
x1 = START_INT+(1-tau)*(END_INT-START_INT);
x2 = START_INT+tau*(END_INT-START_INT);
f_x1 = f(P,x1);
f_x2 = f(P,x2);
while ((abs(END_INT-START_INT)>TOL) && (k2<iter))
if f_x1 < f_x2 %function smaller in x1 than in x2
END_INT = x2; % make interval smaller
x2 = x1; % set new end of interval
x1 = START_INT+(1-tau)*(END_INT-START_INT); % find new beginning
f_x2 = f_x1; % already have value in x1
f_x1 = f(P,x1); % compute new value for new beginning
else
START_INT=x1; % function smaller in x2 than in x1 so set new start indx to x1
x1 = x2; % replace as new start index
x2=START_INT+tau*(END_INT-START_INT); % compute new end index
f_x1= f_x2;
f_x2 = f(P,x2);
end
k2=k2+1;
end
if f_x1 < f_x2
P1 = x1;
epsilon = f_x1;
else
P1 = x2;
epsilon = f_x2;
end
end
EbN0_KaPoissonKnown.m
function data = EbN0_KaPoissonKnown(k, n, epsilon, E_Ka)
% function data = EbN0_KaPoissonKnown(k, n, epsilon, E_Ka)
% Find the minimal required EbN0 (in dB) such that the PUPE is below a
% certain threshold epsilon for a system with the number of active
% users following a Poisson distribution but known.
%
% INPUTS
% k : number of bits per symbol
% n : framelength (number of complex DoFs)
% epsilon : target PUPE
% E_Ka : average number of active users
%
% OUTPUT
% data : store the system parameters and the minimal required EbN0
tic
DEBUG = 1;
%% debugging mode
if DEBUG == 1
k = 100;
n = 15000;
epsilon = 0.1;
E_Ka = 2;
end
%% Poissom PMF of Ka, can be modified to consider other distributions
p_Ka = @(K) poisspdf(K,E_Ka);
%% The range of power to search
EbN0db_lowest = -0.5;
EbN0db_highest = 3;
P_lowest = k*10^(EbN0db_lowest/10)/n;
P_highest = k*10^(EbN0db_highest/10)/n;
%% Function to compute the RCU bound
f_rcu = @(P,P1) RCU_KaRandomKnown(P,P1,k,n,E_Ka,p_Ka);
%% Search for the minimal required EbN0
[eps_RCU, P_RCU] = binary_search(f_rcu, P_lowest, P_highest, epsilon,epsilon/100);
EbN0db_RCU = 10*log10(n*P_RCU/(k));
%% Save the results
sim_time = toc;
data.EbN0db = EbN0db_RCU;
data.E_Ka = E_Ka;
data.eps_est= eps_RCU;
data.epsilon= epsilon;
data.k = k;
data.n = n;
data.sim_time = sim_time;
if DEBUG ~= 1
filename = ['EbN0_KaPoissonKnown_EKa_' num2str(E_Ka) '_epsilon_' num2str(epsilon) '_k_' num2str(k) '_n_' num2str(n) '.mat'];
save(filename, 'data', '-v7.3');
else
keyboard
end
end
RCU_KaPoissonKnown.m
function data = RCU_KaPoissonKnown(k, n, E_Ka, EbN0db)
% function data = RCU_KaPoissonKnown(k, n, E_Ka, EbN0db)
% Compute the RCU bound for the PUPE in Theorem 1 of
%
% [1] Y. Polyanskiy, "A perspective on massive random-access," 2017 IEEE
% International Symposium on Information Theory (ISIT), 2017, pp.
% 2523-2527.
%
% extended to the complex-valued model and the case where Ka follows a
% Poisson distribution and is unknown. See also
%
% [2] K.-H. Ngo, A. Lancho, G. Durisi and A. G. i. Amat, "Unsourced
% Multiple Access With Random User Activity," submitted to IEEE Trans. Inf.
% Theory, Jan. 2022.
%
% INPUTS
% k : number of bits per symbol
% n : framelength (number of complex DoFs)
% E_Ka : average number of active users
% EbN0db : energy per bit in dB
% Note: E_Ka and EbN0db can be vectors.
%
% OUTPUTS
% data : store the system parameters and the computed RCU bound
% P : symbol power budget
% P1 : the actual power used (denoted by P' in [1] and [2])
tic
DEBUG = 0;
%% debugging mode
if DEBUG == 1
k = 128;
n = 19200;
E_Ka = 50;
EbN0db = 0;
end
%% symbol power budget
P_list = k.*10.^(EbN0db./10)./n;
%% initialization
epsilon = zeros(length(EbN0db),length(E_Ka));
P1 = zeros(length(EbN0db),length(E_Ka));
%% Compute the RCU bound
for idxEKa = 1:length(E_Ka)
E_Ka = E_Ka(idxEKa);
p_Ka = @(K) poisspdf(K,E_Ka);
for idxEbN0 = 1:length(EbN0db)
P = P_list(idxEbN0);
f_rcu = @(P,P1) RCU(P,P1,k,n,E_Ka,p_Ka);
[Pe_tmp,P1_tmp] = golden_search(f_rcu,0,P,P/100); % Optimize P'
epsilon(idxEbN0,idxEKa) = Pe_tmp;
P1(idxEbN0,idxEKa) = P1_tmp;
end
end
%% Save the results
sim_time = toc;
data.EbN0db = EbN0db;
data.E_Ka = E_Ka;
data.p_Ka = 'Poisson';
data.epsilon= epsilon;
data.k = k;
data.n = n;
data.P1 = P1;
data.sim_time = sim_time;
if DEBUG ~= 1
filename = ['RCU_KaPoissonKnown_EKa_' num2str(min(E_Ka)) 'to' ...
num2str(max(E_Ka)) '_k_' num2str(k) '_n_' num2str(n) ...
'_EbN0db_' num2str(min(EbN0db)) 'to' num2str(max(EbN0db)) '.mat'];
save(filename, 'data', '-v7.3');
else
keyboard
end
end
RCU_KaFixedKnown.m
function [epsilon, P, P1_opt] = RCU_KaFixedKnown(k,n,Ka,EbN0db)
% function [epsilon] = RCU_KaFixedKnown(k,n,Ka,EbN0db)
% Compute the RCU bound for the PUPE in Theorem 1 of
%
% [1] Y. Polyanskiy, "A perspective on massive random-access," 2017 IEEE
% International Symposium on Information Theory (ISIT), 2017, pp.
% 2523-2527.
%%这个函数是用来计算 Polyanskiy 在他的论文中提到的 RCU（Random Coding Unbounded）界限的
% extended to the complex-valued model.
%
% INPUTS
% k : number of bits per symbol 每符号比特数
% n : framelength (number of complex DoFs) 帧长
% Ka : number of active users 活跃用户数量
% EbN0db : energy per bit in dB 每比特能量（dB）
%
% OUTPUTS
% epsilon : the bound on PUPE given in [1, Eq. (3)]
% epsilon，即ε，表示[1]中给出的PUPE的界线
% P : symbol power budget imposed by the EbN0
% P1_opt : the optimal value of P', the actual power used
%% Example of parameters (for debugging)
% k = 100;
% n = 15000;
% Ka = 100;
% EbN0db = 2;
%% Power and codebook size
P = k*10^(EbN0db/10)/n;
M = 2^k;
%% Initialization
% for the sum over t in Eq. (3)
t_vec = 1:Ka;
% for the optimization over rho and rho1 in E(t)
rho_vec = linspace(0,1,100); % This goes from 0 to 1 in principle
rho1_vec= linspace(0,1,100); % From 0 to 1
rho = repmat(rho_vec,length(rho1_vec),1,length(t_vec));
rho1 = permute(repmat(rho1_vec,length(rho_vec),1,length(t_vec)),[2,1,3]);
t=permute(repmat(t_vec,length(rho1_vec),1,length(rho_vec)),[1,3,2]);
% Some functions for the computation of pt as defined in [1, Th. 1]
R1_f = @(n,M,t) 1./n*log(M)-1./(n*t).*gammaln(t+1);
R2_f = @(n,Ka,t) 1/n*(gammaln(Ka+1)-gammaln(t+1)-gammaln(Ka-t+1));
D_f = @(P1,t,rho,rho1) (P1.*t-1).^2 + 4*P1.*t.*(1+rho.*rho1)./(1+rho);
lambda_f = @(P1,t,rho,rho1) (P1.*t-1+sqrt(D_f(P1,t,rho,rho1)))./(2.*(1+rho1.*rho).*P1.*t); % max(roots_edit([P1.*t.*(rho+1).*(rho.*rho1+1) -P1.*t.*(rho+1) -1]));
mu_f = @(P1,t,rho,rho1) rho.*lambda_f(P1,t,rho,rho1)./(1+P1.*t.*lambda_f(P1,t,rho,rho1));
a_f = @(P1,t,rho,rho1) rho.*log(1+P1.*t.*lambda_f(P1,t,rho,rho1))+log(1+P1.*t.*mu_f(P1,t,rho,rho1));
b_f = @(P1,t,rho,rho1) rho.*lambda_f(P1,t,rho,rho1)- mu_f(P1,t,rho,rho1)./(1+P1.*t.*mu_f(P1,t,rho,rho1));
E0_f = @(P1,t,rho,rho1) rho1.*a_f(P1,t,rho,rho1) + log(1-b_f(P1,t,rho,rho1).*rho1);
Et_f = @(P1,t,rho,rho1,n,M,Ka) squeeze(max(max(-rho.*rho1.*t.*R1_f(n,M,t) - rho1.*R2_f(n,Ka,t) + E0_f(P1,t,rho,rho1))));
pt_f = @(P1,t,rho,rho1,n,M,Ka) exp(-n*Et_f(P1,t,rho,rho1,n,M,Ka));
% number of samples for empirical evaluation of the CDF of It in qt
Niq = 1000;
% for the optimization over P' (here denoted by P1)
P1_vec = linspace(1e-9,P,20);
epsilon = zeros(size(P1_vec));
%% Evaluation of the bound, optimized over P'
% Note: one can perform a golden search for the optimal P'. We do it in the
% evaluation of the RCU bound for Ka random.
for pp = 1:length(P1_vec)
P1 = P1_vec(pp);
% Computation of p0
p0 = nchoosek(Ka,2)/M + Ka*gammainc(n*P/P1,n,'upper');
% Computation of pt
pt = pt_f(P1,t,rho,rho1,n,M,Ka);
% Computation of qt (for t = 1 only, as in [1])
It = zeros(1,1000);
for II = 1:Niq
Zi = sqrt(.5)*(randn(1,n) + 1i*randn(1,n));
codebook = sqrt(.5*P1)*(randn(Ka,n) + 1i*randn(Ka,n));
it = n*log(1+P1) + (sum(abs(repmat(Zi,Ka,1)+codebook).^2,2)./(1+P1)-sum(abs(repmat(Zi,Ka,1)).^2,2));
It(II) = min(it);
end
[prob,gam] = ecdf(It);
qt = min(prob + exp(n*(R1_f(n,M,1)+R2_f(n,Ka,1))-gam));
% Computation of RHS of [1, Eq. (3)] for a given P' < P
epsilon(pp) = t_vec(1)/Ka*min(pt(1),qt) + t_vec(2:end)/Ka*pt(2:end) + p0;
end
% Take the minimum over P'
[epsilon,idx_P1opt] = min([epsilon 1]);
P1_opt = P1_vec(idx_P1opt);
RCU_KaRandomKnown.m
function [epsilon] = RCU_KaRandomKnown(P,P1,k,n,E_Ka,p_Ka)
% function [epsilon] = RCU_KaRandomKnown(P,P1,k,n,E_Ka,p_Ka)
% Compute the RCU bound for the PUPE in Theorem 1 of
%
% [1] Y. Polyanskiy, "A perspective on massive random-access," 2017 IEEE
% International Symposium on Information Theory (ISIT), 2017, pp.
% 2523-2527.
%
% extended to the complex-valued model and the case where Ka is random and
% unknown. See also
%
% [2] K.-H. Ngo, A. Lancho, G. Durisi and A. G. i. Amat, "Unsourced
% Multiple Access With Random User Activity," submitted to IEEE Trans. Inf.
% Theory, Jan. 2022.
%
% INPUTS
% P : symbol power budget
% P1 : the actual power used (denoted by P' in [1] and [2])
% k : number of bits per symbol
% n : framelength (number of complex DoFs)
% E_Ka : average number of active users
% p_Ka : PMF of the number of active users Ka
%
% OUTPUTS
% epsilon : (extended) bound on PUPE given in [1, Eq. (3)]
% codebook size
M = 2^k;
% number of samples for empirical evaluation of the CDF of It in qt
Niq = 1000;
%% Computation of p0 (with some additional terms w.r.t. p0 in [1] due to the randomness of Ka. See [2].)
% find K_l and K_u such that the probability that Ka is not in [K_l:K_u] is small
K_l = floor(E_Ka); K_u = ceil(E_Ka);
while p_Ka(K_l-1) > 1e-9
K_l = K_l - 1;
end
K_l = max(K_l,0);
while p_Ka(K_u+1) > 1e-9
K_u = K_u + 1;
end
p01 = 0;
for Ka_tmp = K_l:K_u
p01_tmp = 1;
for ii = 1:Ka_tmp-1
p01_tmp = p01_tmp*(1-ii/M);
end
p01 = p01 + p_Ka(Ka_tmp)*p01_tmp;
% p01 = p01 + p_Ka(Ka_tmp)*nchoosek(Ka_tmp,2)/M;
end
p0 = 2 - p01 - sum(p_Ka(K_l:K_u)) + E_Ka*gammainc(n*P/P1,n,'upper');
%% Initialize epsilon to p0
epsilon = p0;
%% Compute the RCU bound, averaged over the distribution of Ka
parfor Ka = max(K_l,1):K_u
%% Computation of pt
% Initialize for the sum over t and the optimization over rho and rho1
rho_vec = linspace(0,1,100);
rho1_vec= linspace(0,1,100);
t_vec = 1:Ka;
rho = repmat(rho_vec,length(rho1_vec),1,length(t_vec));
rho1 = permute(repmat(rho1_vec,length(rho_vec),1,length(t_vec)),[2,1,3]);
t = permute(repmat(t_vec,length(rho1_vec),1,length(rho_vec)),[1 3 2]);
% Some functions for the computation of pt as defined in [1, Th. 1]
R1_f = @(n,M,t) 1./n*log(M)-1./(n*t).*gammaln(t+1);
R2_f = @(n,Ka,t) 1/n*(gammaln(Ka+1)-gammaln(t+1)-gammaln(Ka-t+1));
D_f = @(P1,t,rho,rho1) (P1.*t-1).^2 + 4*P1.*t.*(1+rho.*rho1)./(1+rho);
lambda_f = @(P1,t,rho,rho1) (P1.*t-1+sqrt(D_f(P1,t,rho,rho1)))./(2.*(1+rho1.*rho).*P1.*t); % max(roots_edit([P1.*t.*(rho+1).*(rho.*rho1+1) -P1.*t.*(rho+1) -1]));
mu_f = @(P1,t,rho,rho1) rho.*lambda_f(P1,t,rho,rho1)./(1+P1.*t.*lambda_f(P1,t,rho,rho1));
a_f = @(P1,t,rho,rho1) rho.*log(1+P1.*t.*lambda_f(P1,t,rho,rho1))+log(1+P1.*t.*mu_f(P1,t,rho,rho1));
b_f = @(P1,t,rho,rho1) rho.*lambda_f(P1,t,rho,rho1)- mu_f(P1,t,rho,rho1)./(1+P1.*t.*mu_f(P1,t,rho,rho1));
E0_f = @(P1,t,rho,rho1) rho1.*a_f(P1,t,rho,rho1) + log(1-b_f(P1,t,rho,rho1).*rho1);
Et_f = @(P1,t,rho,rho1,n,M,Ka) squeeze(max(max(-rho.*rho1.*t.*R1_f(n,M,t) - rho1.*R2_f(n,Ka,t) + E0_f(P1,t,rho,rho1))));
pt_f = @(P1,t,rho,rho1,n,M,Ka) exp(-n*Et_f(P1,t,rho,rho1,n,M,Ka));
% Compute pt
pt = pt_f(P1,t,rho,rho1,n,M,Ka);
%% Computation of qt (for t = 1 only, as in [1])
qt = 1;
% compute qt for Ka < 50 only due to complexity issue
if Ka < 50
It = zeros(1,Niq);
for II = 1:Niq
Zi = sqrt(.5)*(randn(1,n) + 1i*randn(1,n));
codebook = sqrt(.5*P1)*(randn(Ka,n) + 1i*randn(Ka,n));
it = n*log(1+P1) + (sum(abs(repmat(Zi,Ka,1)+codebook).^2,2)./(1+P1)-sum(abs(repmat(Zi,Ka,1)).^2,2));
It(II) = min(it);
end
[prob,gam] = ecdf(It);
qt = min(prob + exp(n*(R1_f(n,M,1)+R2_f(n,Ka,1))-gam));
end
%% Take the min of pt and qt
pt(1) = min(pt(1),qt);
%% Compute epsilon, the RHS of [1, Eq. (3)]
epsilon = epsilon + (t_vec/Ka*pt)*feval(p_Ka,Ka);
end
end
RCU_KaUnknown
golden_search_P1_MDFA.m
function [eps_MD, eps_FA, P1] = golden_search_P1_MDFA(f, START_INT, END_INT, TOL, type, weight_MD, weight_FA)
% Minimize a weighted sum of eps_MD and eps_FA over P1 \in [0,P], where
% eps_MD and eps_FA are given by f(P,P1), which is the RCU bound in Theorem
% 1 of
%
% [1] K.-H. Ngo, A. Lancho, G. Durisi and A. G. i. Amat, "Unsourced
% Multiple Access With Random User Activity," submitted to IEEE Trans. Inf.
% Theory, Jan. 2022.
if strcmp(type,'max')
g = @(fMD,fFA) max(fMD,fFA);
elseif strcmp(type,'weighted')
g = @(fMD,fFA) fMD*weight_MD + fFA*weight_FA;
end
P = END_INT; % power budget
tau = double((sqrt(5)-1)/2); % golden proportion coefficient, around 0.618
iter = 20; % maximum number of iterations
k2 = 0; % iteration index
x1 = START_INT+(1-tau)*(END_INT-START_INT);
x2 = START_INT+tau*(END_INT-START_INT);
[fMD_x1,fFA_x1] = f(P,x1);
[fMD_x2,fFA_x2] = f(P,x2);
while ((abs(END_INT-START_INT)>TOL) && (k2<iter))
if g(fMD_x1, fFA_x1) < g(fMD_x2, fFA_x2) %function smaller in x1 than in x2
END_INT = x2; %make interval smaller
x2 = x1; %set new end of interval
x1 = START_INT+(1-tau)*(END_INT-START_INT); %find new beginning
fMD_x2 = fMD_x1;%already have value in x1
fFA_x2 = fFA_x1;
[fMD_x1,fFA_x1] = f(P,x1); %%compute new value for new beginning
else
START_INT=x1; %function smaller in x2 than in x1 so set new start indx to x1
x1 = x2; %replace as new start index
x2=START_INT+tau*(END_INT-START_INT); %compute new end index
fMD_x1= fMD_x2;
fFA_x1 = fFA_x2;
[fMD_x2,fFA_x2] = f(P,x2);
end
k2=k2+1;
end
% Return the solution
if g(fMD_x1, fFA_x1) < g(fMD_x2, fFA_x2)
P1 = x1;
eps_MD = fMD_x1;
eps_FA = fFA_x1;
else
P1 = x2;
eps_MD = fMD_x2;
eps_FA = fFA_x2;
end
end
RCU_floor_KaRandomUnknown.m
function [floor_MD,floor_FA] = RCU_floor_KaRandomUnknown(rad_l,rad_u,k,n,E_Ka,p_Ka)
% function [floor_MD,floor_FA] = RCU_floor(rad_l,rad_u,k,n,E_Ka,p_Ka)
% Compute the error floors for the RCU bounds of the misdetection and false
% alarm probabilities, characterized in Theorem 3 of
%
% [1] K.-H. Ngo, A. Lancho, G. Durisi and A. G. i. Amat, "Unsourced
% Multiple Access With Random User Activity," submitted to IEEE Trans. Inf.
% Theory, Jan. 2022.
%
% for a system with random and unknown number of active users.
%
% INPUTS
% rad_l : lower decoding radius
% rad_u : upper decoding radius
% k : number of bits per symbol
% n : framelength (number of complex DoFs)
% E_Ka : average number of active users
% p_Ka : PMF of the number of active users Ka
%
% OUTPUTS
% floor_MD, floor_FA : the error floors
% codebook size
M = 2^k;
%% Computation of \bar{p}
K_l = floor(E_Ka); K_u = ceil(E_Ka);
while p_Ka(K_l-1) > 1e-12
K_l = K_l - 1;
end
K_l = max(K_l,0);
while p_Ka(K_u+1) > 1e-12
K_u = K_u + 1;
end
p01 = 0;
for Ka_tmp = K_l:K_u
p01_tmp = 1;
for ii = 1:Ka_tmp-1
p01_tmp = p01_tmp*(1-ii/M);
end
p01 = p01 + p_Ka(Ka_tmp)*p01_tmp;
end
pbar = 1 - p01 + 1 - sum(p_Ka(K_l:K_u));
%% Initialize floor_MD and floor_FA to be \bar{p}
floor_MD = pbar;
floor_FA = pbar;
%% The sum over Ka
parfor Ka = K_l:K_u
% Compute \xi(Ka,Ka')
KaEst_thres = @(Ka,KaEst) n.*log(Ka./KaEst)./(Ka./KaEst-1);
P_Ka_KaEst_list_a = gammainc(KaEst_thres(Ka,K_l:Ka-1),n,'lower');
P_Ka_KaEst_list_b = gammainc(KaEst_thres(Ka,Ka+1:K_u),n,'upper');
P_Ka_Ka = 1 - max([P_Ka_KaEst_list_a P_Ka_KaEst_list_b]);
P_Ka_KaEst_list = [P_Ka_KaEst_list_a'; P_Ka_Ka; P_Ka_KaEst_list_b'];
%% The sum over Ka'
for KaEst = K_l:K_u
if KaEst == 0
P_Ka_KaEst = 0;
else
P_Ka_KaEst = P_Ka_KaEst_list(KaEst - K_l + 1);
end
KaEst_l = max(KaEst - rad_l,K_l);
KaEst_u = min(KaEst + rad_u,K_u);
if Ka > 0
floor_MD = floor_MD + feval(p_Ka, Ka)*max(Ka-KaEst_u,0)*P_Ka_KaEst/Ka;
end
Mrx = (Ka + max(KaEst_l-Ka,0) - max(Ka-KaEst_u,0)); % number of decoded codewords
if Mrx > 0
floor_FA = floor_FA + feval(p_Ka, Ka)*max(KaEst_l-Ka,0)*P_Ka_KaEst/Mrx;
end
end
end
floor_MD = min(floor_MD,1);
floor_FA = min(floor_FA,1);
end
binary_search_P_MDFA.m
function [eps_MD, eps_FA, P,P1] = binary_search_P_MDFA(f,x1,x2,TARGET_MD,TARGET_FA,TOL,fixP1)
% Search the value of P such that eps_MD \in [TARGET_MD - TOL, TARGET_MD]
% and eps_FA \in [TARGET_FA - TOL, TARGET_FA] where eps_MD and eps_FA are
% obtained from a f(P,P1), which is the RCU bound in Theorem 1 of
%
% [1] K.-H. Ngo, A. Lancho, G. Durisi and A. G. i. Amat, "Unsourced
% Multiple Access With Random User Activity," submitted to IEEE Trans. Inf.
% Theory, Jan. 2022.
iter= 20; % maximum number of iterations
k1=0; % number of iterations
%% We can either optimize over P1 or choose a fixed value of P1. The latter
%% is faster. Our experiments suggest that the optimal P1 is around 0.96*P
fracP1 = 0.96;
if fixP1 == 0
if TARGET_MD == TARGET_FA
search_P1 = @(x) golden_search_P1_MDFA(f,1e-9,x,x/100,'max',[],[]);
else
% Calculate the weights based on the target MD and FA probabilities
weight_MD = 1/(1 + TARGET_MD/TARGET_FA);
weight_FA = 1/(1 + TARGET_FA/TARGET_MD);
search_P1 = @(x) golden_search_P1_MDFA(f,1e-9,x,x/100,'weighted',weight_MD, weight_FA);
end
else
search_P1 = @(x) compute_MDFA_fixedP1(f, x, fracP1);
end
%% Search for P
[eps_MD_tmp,eps_FA_tmp,P1_tmp] = search_P1(x2);
[eps_MD_tmp1,eps_FA_tmp1,P1_tmp1] = search_P1(x1);
if x1 == x2
P = x1;
eps_MD = eps_MD_tmp;
eps_FA = eps_FA_tmp;
P1 = P1_tmp;
elseif eps_MD_tmp > TARGET_MD || eps_FA_tmp > TARGET_FA
warning('Impossible to achieve the target within the given range of power :( ');
P = x2;
eps_MD = eps_MD_tmp;
eps_FA = eps_FA_tmp;
P1 = P1_tmp;
elseif eps_MD_tmp1 < TARGET_MD || eps_FA_tmp1 < TARGET_FA
warning('All power values in the range can achieve the target :) ');
P = x1;
eps_MD = eps_MD_tmp1;
eps_FA = eps_FA_tmp1;
P1 = P1_tmp1;
else
[fx_MD,fx_FA,P1_tmp] = search_P1((x1+x2)/2);
while ~((TARGET_MD >= fx_MD && fx_MD >= TARGET_MD - TOL && TARGET_FA >= fx_FA) || ...
(TARGET_FA >= fx_FA && fx_FA >= TARGET_FA - TOL && TARGET_MD >= fx_MD) ...
|| (k1>iter))
if k1 > 0
[fx_MD,fx_FA,P1_tmp] = search_P1((x1+x2)/2);
end
if TARGET_MD > fx_MD && TARGET_FA > fx_FA
x2 = (x1+x2)/2; % set new end of interval
else
x1 = (x1+x2)/2; % replace as new start index
end
k1 = k1+1;
end
eps_MD = fx_MD;
eps_FA = fx_FA;
P = x2;
P1 = P1_tmp;
end
end
function [eps_MD, eps_FA, P1] = compute_MDFA_fixedP1(f, P, fracP1)
P1 = P*fracP1;
[eps_MD, eps_FA] = f(P,P1);
end
EbN0_KaPoissonUnknown.m
function data = EbN0_KaPoissonUnknown(k, n, epsilon_MD, epsilon_FA, E_Ka, rad_l, rad_u)
% function data = EbN0_KaPoissonUnknown(k, n, epsilon_MD, epsilon_FA, E_Ka)
% Find the minimal required EbN0 (in dB) such that the misdetection and
% false alarm probabilities are below the thresholds epsilon_MD and
% epsilon_FA, respectively, for a system with the number of active
% users following a Poisson distribution and unknown.
%
% INPUTS
% k : number of bits per symbol
% n : framelength (number of complex DoFs)
% epsilon_MD : target misdetection probability
% epsilon_FA : target false alarm probability
% E_Ka : average number of active users
% rad_l : lower decoding radius
% rad_u : upper decoding radius
%
% OUTPUT
% data : store the system parameters and the minimal required EbN0
tic
DEBUG = 0;
%% debugging mode
if DEBUG == 1
k = 128; % Number of bits
n = 19200; % Frame length
epsilon_MD = .1; % Per-user error probability
epsilon_FA = .1; % Per-user error probability
E_Ka = 50;
rad_l = 0;
rad_u = 0;
end
%% Poissom PMF of Ka, can be modified to consider other distributions
p_Ka = @(K) poisspdf(K,E_Ka);
%% The range of power to search
EbN0db_lowest = -2;
EbN0db_highest = 15;
P_lowest = k*10^(EbN0db_lowest/10)/n;
P_highest = k*10^(EbN0db_highest/10)/n;
%% The below can be used to find suitable decoding radii
% rad_l = 0; rad_u = 0;
%
% [floor_MD,floor_FA] = RCU_floor(rad_l,rad_u,k,n,E_Ka,p_Ka);
% while floor_MD > epsilon_MD || floor_FA > epsilon_FA
% if floor_MD > epsilon_MD
% rad_u = rad_u + 1;
% end
% if floor_FA > epsilon_FA
% rad_l = rad_l + 1;
% end
% [floor_MD,floor_FA] = RCU_floor(rad_l,rad_u,k,n,E_Ka,p_Ka);
% end
% if (epsilon_MD < 1e-1) || (epsilon_FA < 1e-1)
% rad_u = rad_u + 2;
% rad_l = rad_l + 2;
% end
%% Function to compute the RCU bound
f_rcu = @(P,P1) RCU_KaRandomUnknown(P,P1,rad_l,rad_u,k,n,E_Ka,p_Ka);
%% Search for the minimal required EbN0
[eps_RCU_MD, eps_RCU_FA, P_RCU,P1] = binary_search_P_MDFA(f_rcu, P_lowest, P_highest,...
epsilon_MD,epsilon_FA,min(epsilon_MD,epsilon_FA)/100,0);
EbN0db_RCU = 10*log10(n*P_RCU/k);
%% Save the results
sim_time = toc;
data.EbN0db = EbN0db_RCU;
data.E_Ka = E_Ka;
data.p_Ka = 'Poisson';
data.eps_est_MD = eps_RCU_MD;
data.eps_est_FA = eps_RCU_FA;
data.epsilon_MD = epsilon_MD;
data.epsilon_FA = epsilon_FA;
data.rad_l = rad_l;
data.rad_u = rad_u;
data.k = k;
data.n = n;
data.P1 = P1;
data.sim_time = sim_time;
if DEBUG ~= 1
filename = ['EbN0_KaPoissonUnknown_EKa_' num2str(E_Ka) '_epsilonMD_' ...
num2str(epsilon_MD) '_epsilonFA_' num2str(epsilon_FA) ...
'_k_' num2str(k) '_n_' num2str(n) '.mat'];
save(filename, 'data', '-v7.3');
else
keyboard
end
end
RCU_KaPoissonUnknown.m
function data = RCU_KaPoissonUnknown(k, n, E_Ka_list, EbN0db, rad_l_list, rad_u_list)
% function data = RCU_KaPoissonUnknown(k, n, E_Ka_list, EbN0db, rad_l_list, rad_u_list)
% Compute the RCU bound for the PUPE for a system with number of users
% unknown and following a Poission distribution. The bound is given in
% Theorem 1 of
%
% [1] K.-H. Ngo, A. Lancho, G. Durisi and A. G. i. Amat, "Unsourced
% Multiple Access With Random User Activity," submitted to IEEE Trans. Inf.
% Theory, Jan. 2022.
%
% INPUTS
% k : number of bits per symbol
% n : framelength (number of complex DoFs)
% E_Ka_list : set of values for the average number of active users
% EbN0db : energy per bit in dB
% rad_l_list : set of values for the lower decoding radius
% rad_u_list : set of corresponding upper decoding radius
% Note: rad_l_list and rad_u_list are of the same length
%
% OUTPUTS
% data : store the system parameters and the computed RCU bound
tic
DEBUG = 0;
%% debugging mode
if DEBUG == 1
k = 128;
n = 19200;
EbN0db = 4;
E_Ka_list = 50;
rad_l_list = [0];
rad_u_list = [0];
end
%% codebook size
M = 2^k;
%% symbol power budget
P_list = k.*10.^(EbN0db./10)./n;
%% initialization
p_MD = ones(length(EbN0db),length(E_Ka_list),length(rad_l_list));
p_FA = ones(length(EbN0db),length(E_Ka_list),length(rad_l_list));
P1 = zeros(length(EbN0db),length(E_Ka_list),length(rad_l_list));
%% Compute the RCU
for idxEKa = 1:length(E_Ka_list)
E_Ka = E_Ka_list(idxEKa);
p_Ka = @(K) poisspdf(K,E_Ka);
for idxRad = 1:length(rad_l_list)
rad_l = rad_l_list(idxRad); % lower decoding radius
rad_u = rad_u_list(idxRad); % upper decoding radius
for idxEbN0 = 1:length(EbN0db)
P = P_list(idxEbN0);
f_rcu = @(P,P1) RCU_KaRandomUnknown(P,P1,rad_l,rad_u,k,n,E_Ka,p_Ka);
% Optimize P1 to minimize max(p_MD,p_FA)
[p_MD_tmp,p_FA_tmp,P1_tmp] = golden_search_P1_MDFA(f_rcu,0,P,P/200,'max',[],[]);
% To minimize a weighted sum of p_MD and p_FA, use the below
% weight_MD = 1;
% weight_FA = 1;
% [p_MD_tmp,p_FA_tmp,P1_tmp] = golden_search_P1_MDFA(f_rcu,0,P,P/100,'weighted',weight_MD,weight_FA);
p_MD(idxEbN0,idxEKa,idxRad) = p_MD_tmp;
p_FA(idxEbN0,idxEKa,idxRad) = p_FA_tmp;
P1(idxEbN0,idxEKa,idxRad) = P1_tmp;
end
end
end
p_MD = squeeze(p_MD);
p_FA = squeeze(p_FA);
P1 = squeeze(P1);
%% Save the results
sim_time = toc;
data.EbN0db = EbN0db;
data.E_Ka = E_Ka_list;
data.p_Ka = 'Poisson';
data.k = k;
data.n = n;
data.rad_lower = rad_l;
data.rad_upper = rad_u;
data.p_MD = p_MD;
data.p_FA = p_FA;
data.P1 = P1;
data.sim_time = sim_time;
if DEBUG ~= 1
filename = ['RCU_KaPoissonUnknown_EKa_' num2str(min(E_Ka_list)) 'to' ...
num2str(max(E_Ka_list)) '_k_' num2str(k) '_n_' num2str(n) ...
'_radL_' sprintf('%d', rad_l_list) ...
'_radU_' sprintf('%d', rad_u_list) ...
'_EbN0db_' num2str(min(EbN0db)) 'to' num2str(max(EbN0db)) '.mat'];
save(filename, 'data', '-v7.3');
else
keyboard
end
end
RCU_KaRandomUnknown.m
function [eps_MD,eps_FA] = RCU_KaRandomUnknown(P,P1,rad_l,rad_u,k,n,E_Ka,p_Ka)
% function [eps_MD,eps_FA] = RCU_KaRandomUnknown(P,P1,rad_l,rad_u,k,n,E_Ka,p_Ka)
% Compute the RCU bound for the PUPE in Theorem 1 of
%
% [1] K.-H. Ngo, A. Lancho, G. Durisi and A. G. i. Amat, "Unsourced
% Multiple Access With Random User Activity," submitted to IEEE Trans. Inf.
% Theory, Jan. 2022.
%
% for a system with random and unknown number of active users.
%
% INPUTS
% P : symbol power budget
% P1 : the actual power used (denoted by P' in [1] and [2])
% rad_l : lower decoding radius
% rad_u : upper decoding radius
% k : number of bits per symbol
% n : framelength (number of complex DoFs)
% E_Ka : average number of active users
% p_Ka : PMF of the number of active users Ka
%
% OUTPUTS
% eps_MD, eps_FA : RCU bounds on the misdetection and false alarm
% probabilities, respectively
% codebook size
M = 2^k;
% number of samples for empirical evaluation of the CDF of It in qt
Niq = 1000;
%% Computation of \tilde{p}
% Compute the error floors [1, Th. 3]
[floor_MD,floor_FA] = RCU_floor_KaRandomUnknown(rad_l,rad_u,k,n,E_Ka,p_Ka);
% Find K_l and K_u such that Pr[Ka \notin [K_l:K_u]] is small
K_l = floor(E_Ka); K_u = ceil(E_Ka);
while p_Ka(K_l-1) > max(.0001*min(floor_MD,floor_FA),1e-9)
K_l = K_l - 1;
end
K_l = max(K_l,0);
while p_Ka(K_u+1) > max(.0001*min(floor_MD,floor_FA),1e-9)
K_u = K_u + 1;
end
p01 = 0;
for Ka_tmp = K_l:K_u
p01_tmp = 1;
for ii = 1:Ka_tmp-1
p01_tmp = p01_tmp*(1-ii/M);
end
p01 = p01 + p_Ka(Ka_tmp)*p01_tmp;
end
ptilde = 1 - p01 + 1 - sum(p_Ka(K_l:K_u)) + E_Ka*gammainc(n*P/P1,n,'upper');
if ptilde >= 1
eps_MD = 1;
eps_FA = 1;
return
end
%% Initialize eps_MD and eps_FA to be p0
eps_MD = ptilde;
eps_FA = ptilde;
%% The expectation w.r.t. Ka
parfor Ka = K_l:K_u
%% Compute the term \xi(Ka,Ka') using [1, Th. 2]
% The function \zeta for ML esimator of Ka
KaEst_thres = @(Ka,KaEst,P1) n.*(log(1+Ka.*P1) - log(1+KaEst.*P1))./((1+Ka.*P1)./(1+KaEst.*P1)-1);
% For energy-based estimator of Ka, use the below
% KaEst_thres = @(Ka,KaEst,P1) n./(1+Ka.*P1).*(1+(Ka+KaEst).*P1/2);
% Option 1: compute the exact term defined in [1]
% tic
% P_Ka_KaEst_list = zeros(Ka_u - Ka_l + 1,1);
% for idx = 1:Ka_u - Ka_l + 1
% KaEst = Ka_l + idx - 1;
% P_Ka_KaEst_list_a = gammainc(KaEst_thres(Ka,Ka_l:KaEst-1,P1),n,'lower');
% P_Ka_KaEst_list_b = gammainc(KaEst_thres(Ka,KaEst+1:Ka_u,P1),n,'upper');
% P_Ka_KaEst_list(idx) = 1 - max([P_Ka_KaEst_list_a P_Ka_KaEst_list_b]);
% end
% toc
% Option 2: slight relaxation that is faster
% tic
P_Ka_KaEst_list_a = gammainc(KaEst_thres(Ka,K_l:Ka-1,P1),n,'lower');
P_Ka_KaEst_list_b = gammainc(KaEst_thres(Ka,Ka+1:K_u,P1),n,'upper');
P_Ka_Ka = 1 - max([P_Ka_KaEst_list_a P_Ka_KaEst_list_b]);
P_Ka_KaEst_list = [P_Ka_KaEst_list_a'; P_Ka_Ka; P_Ka_KaEst_list_b'];
% toc
%% The expectation w.r.t. Ka'
for KaEst = K_l:K_u
% \xi(Ka,Ka')
xi_Ka_KaEst = P_Ka_KaEst_list(KaEst - K_l + 1);
% \underbar{Ka'}
KaEst_l = max(KaEst - rad_l,K_l);
% \overbar{Ka'}
KaEst_u = min(KaEst + rad_u,K_u);
% Sets of possible values of t and t'
t_vec = 0:min(min(KaEst_u,Ka),M-KaEst_l-max(Ka-KaEst_u,0));
t1_l = max(Ka-KaEst_u,0) - max(Ka-KaEst_l,0);
t1_u = max(KaEst_u-Ka,0) - max(KaEst_l-Ka,0);
t1_vec = zeros(length(t_vec),KaEst_u-KaEst_l+1);
for idxT = 1:length(t_vec)
t1_vec(idxT,:) = t1_l+t_vec(idxT):t1_u+t_vec(idxT);
end
%% Computation of pt:
% Definitions of functions R_1 and R_2
if M-max(Ka,KaEst_l) <= 2^20
R1_f = @(n,Ka,K1,M,t1) (gammaln(M-max(Ka,K1) + 1) - gammaln(M-max(Ka,K1)-t1+ 1) - gammaln(t1+1))/n;
else
R1_f = @(n,Ka,K1,M,t1) t1./n.*log(M-max(Ka,K1))-1./n.*gammaln(t1+1);
end
R2_f = @(n,Ka,K1,t) 1/n*(gammaln(min(Ka,K1)+1)-gammaln(t+1)-gammaln(min(Ka,K1)-t+1));
% First, consider t1 > 0
t1_vec2 = t1_vec;
t1_vec2(t1_vec2<0) = 0;
% Initialize for the optimization over rho and rho1
rho_vec = linspace(1e-9,1,50);
rho1_vec = linspace(1e-9,1,50);
rho = permute(repmat(rho_vec,length(rho1_vec),1,length(t_vec),size(t1_vec,2)),[2,1,3,4]);
rho1 = repmat(rho1_vec,length(rho_vec),1,length(t_vec),size(t1_vec,2));
t = permute(repmat(t_vec,length(rho_vec),1,length(rho1_vec),size(t1_vec,2)),[1,3,2,4]);
t1 = permute(repmat(t1_vec2,1,1,length(rho_vec),length(rho1_vec)),[3,4,1,2]);
% Precompute some quantities
P2 = 1+(max(Ka-KaEst_u,0)+max(KaEst_l-Ka,0)).*P1;
P3 = P1.*(t1+rho.*t);
P4 = rho.*rho1.*P3.*P2;
P5 = rho1.*(rho-1).*t1.*P1;
% Find optimal lambda that maximizes E0(t,t')
c1_f = - P3.*P4.*P5 - (rho1+1).*t1.*P1.*P3.*P4;
c2_f = P5.*(P3.^2 - P4) + rho1.*(t1.*P1.*P3.^2 - P3.*P4) ...
- (2.*t1.*P1 + P3).*P4;
c3_f = 2.*P3.*P5 + rho1.*(P3.^2 + t1.*P1.*P3) - 2.*P4;
c4_f = P5 + rho1.*P3;
Delta0_f = c2_f.^2 - 3.*c1_f.*c3_f;
Delta1_f = 2.*c2_f.^3 - 9.*c1_f.*c2_f.*c3_f + 27.*c1_f.^2.*c4_f;
Q_f = ((Delta1_f + sqrt(Delta1_f.^2 - 4.*Delta0_f.^3))./2).^(1/3);
lambda = real(-(c2_f + Q_f + Delta0_f./Q_f)./3./c1_f);
% Compute E0(t,t')
E0_f = rho1.*(rho-1).*log(1+lambda.*P1.*t1) ...
+ (rho1-1).*log(1+lambda.*P3) ...
+ log(1+lambda.*P3 - lambda.^2.*P4);
% Compute E(t,t')
Ett_f = reshape(max(max(-rho.*rho1.*R1_f(n,Ka,KaEst_l,M,t1) ...
- rho1.*R2_f(n,Ka,KaEst_u,t) + E0_f)),[length(t_vec) size(t1_vec,2)]);
% Compute p_{t,t'}
ptt = min(exp(-n.*Ett_f),1);
% Now, consider the case t1 = 0, where the optimal lambda is the root
% of a quadratic function
rho = permute(repmat(rho_vec,length(rho1_vec),1,length(t_vec)),[2,1,3]);
rho1 = repmat(rho1_vec,length(rho_vec),1,length(t_vec));
t = permute(repmat(t_vec,length(rho_vec),1,length(rho1_vec)),[1,3,2]);
P3 = P1.*rho.*t;
P4 = rho.*rho1.*P3.*P2;
c2_f = -(rho1+1).*P3.*P4;
c3_f = rho1.*P3.^2 - 2.*P4;
c4_f = rho1.*P3;
Delta = c3_f.^2 - 4.*c2_f.*c4_f;
lambda = -(c3_f+sqrt(Delta))./2./c2_f;
E0_f = (rho1-1).*log(1+lambda.*P3) + log(1+lambda.*P3 - lambda.^2.*P4);
Ett_f = squeeze(max(max(- rho1.*R2_f(n,Ka,KaEst_u,t) + E0_f)));
ptt_t1zero = min(exp(-n.*Ett_f),1);
% Combine the cases t1 > 0 and t1 = 0
ptt(t1_vec > min(KaEst_u-max(KaEst_l-Ka,0),M-max(KaEst_l,Ka))) = 0;
ptt(t1_vec < 0) = 0;
for idxT = 1:length(t_vec)
for idxT1 = 1:size(t1_vec,2)
if t1_vec(idxT,idxT1) == 0
ptt(idxT,idxT1) = ptt_t1zero(idxT);
end
end
end
% Compute pt
pt = min(sum(ptt,2),1);
%% Computation of qt for t = 1:
qt = 1;
qtt = 1;
if t_vec(end) >= 1
t1_set = t1_vec2(2,:);
% Compute qt for t = 1 and Ka <= 50 only due to complexity issue
if t_vec(1) <= 1 && Ka <= 150
It = zeros(1,Niq);
for II = 1:Niq
Zi = sqrt(.5)*(randn(1,n) + 1i*randn(1,n));
DefaultFA = sqrt(.5*max(KaEst_l-Ka,0)*P1)*(randn(1,n) + 1i*randn(1,n));
codebook = sqrt(.5*P1)*(randn(min(Ka,KaEst_u),n) + 1i*randn(min(Ka,KaEst_u),n));
it = n*log(1+(1+max(Ka-KaEst_u,0))*P1) + ...
(sum(abs(repmat(Zi+DefaultFA,min(Ka,KaEst_u),1)+codebook).^2,2)./...
(1+(1+max(Ka-KaEst_u,0))*P1)-sum(abs(repmat(Zi,min(Ka,KaEst_u),1)).^2,2));
It(II) = min(it);
end
[prob,gam] = ecdf(It);
qt = min(prob + sum(exp(repmat(n*(R1_f(n,Ka,KaEst_l,M,t1_set)+R2_f(n,Ka,KaEst_u,1)),length(gam),1)-gam),2));
qtt = min(repmat(prob,1,length(t1_set)) + exp(n*(R1_f(n,Ka,KaEst_l,M,t1_set)+R2_f(n,Ka,KaEst_u,1))-gam));
end
% Take the min of pt and qt
pt(2) = min(pt(2),qt);
ptt(2,:) = min(ptt(2,:),qtt);
end
% Take the min of pt, qt, and xi(Ka,Ka')
pt = min(pt,xi_Ka_KaEst);
ptt = min(ptt,xi_Ka_KaEst);
%% Computation epsilon_MD and epsilon_FA for a given P'< P
if Ka > 0
eps_MD = eps_MD + feval(p_Ka, Ka)*(sum((t_vec + max(Ka-KaEst_u,0))*pt)/Ka);
end
t_vec_2 = repmat(t_vec,size(t1_vec,2),1).';
Mrx = (Ka-t_vec_2+t1_vec + max(KaEst_l-Ka,0) - max(Ka-KaEst_u,0)); % number of decoded codewords
Mrx(Mrx == 0) = inf;
eps_FA = eps_FA + feval(p_Ka, Ka)*sum(sum(ptt.*(t1_vec + max(KaEst_l-Ka,0))./Mrx,1));
end
% if eps_MD >= 1 && eps_FA >= 1
% break;
% end
end
eps_MD = min(eps_MD,1);
eps_FA = min(eps_FA,1);
end
SA-MPR
EbN0_SAMPR_KaPoissonKnown.m
function data = EbN0_SAMPR_KaPoissonKnown(k, n, epsilon, E_Ka, SlotIdxCoding)
% function data = EbN0_SAMPR_KaPoissonKnown(k, n, epsilon, E_Ka, SlotIdxCoding)
% Find the minimal required EbN0 (in dB) such that the PUPE achieved with
% slotted ALOHA (SA) with multi-packet reception (MPR) is below a
% certain threshold epsilon for a system with the number of active
% users following a Poisson distribution, but known.
%
% INPUTS
% k : number of bits per symbol
% n : framelength (number of complex DoFs)
% epsilon : target PUPE
% E_Ka : average number of active users
% SlotIdxCoding : 1 if slot-index coding is employed, 0 otherwise
%
% OUTPUT
% data : store the system parameters and the minimal required EbN0
addpath ../RCU_KaKnown
tic
DEBUG = 0;
%% debugging mode
if DEBUG == 1
k = 128;
n = 19200;
epsilon = .1;
E_Ka = 50;
end
%% Poissom PMF of Ka, can be modified to consider other distributions
p_Ka = @(K) poisspdf(K,E_Ka);
%% The range of power to search
EbN0db_lowest = -1;
EbN0db_highest = 15;
P_lowest = k*10^(EbN0db_lowest/10)/n;
P_highest = k*10^(EbN0db_highest/10)/n;
%% Function to compute the RCU bound
if SlotIdxCoding == 1
f_rcu = @(P,P1,L) RCU_KaRandomKnown(P*L,P1*L,k-floor(log2(L)),...
ceil(n/L),E_Ka/L,@(Ksa) PMF_Ksa(p_Ka,E_Ka,L,Ksa));
else
f_rcu = @(P,P1,L) RCU_KaRandomKnown(P*L,P1*L,k,ceil(n/L),E_Ka/L,...
@(Ksa) PMF_Ksa(p_Ka,E_Ka,L,Ksa));
end
%% Search for the minimal required EbN0
% The RCU is optimized over P1 and L
[eps_RCU,P_RCU,P1,L] = binary_search_P_SAMPR_KaKnown(f_rcu, P_lowest, P_highest,...
epsilon,epsilon/100,E_Ka);
EbN0db_RCU = 10*log10(n*P_RCU/k);
%% Save the results
sim_time = toc;
data.EbN0db = EbN0db_RCU;
data.E_Ka = E_Ka;
data.p_Ka = 'Poisson';
data.eps_RCU = eps_RCU;
data.epsilon = epsilon;
data.k = k;
data.n = n;
data.P1 = P1;
data.nSlot = L;
data.SlotIdxCoding = SlotIdxCoding;
data.sim_time = sim_time;
if DEBUG ~= 1
if SlotIdxCoding == 1
filename = ['EbN0_SAMPR_SlotIdxCoding_KaPoissonKnown_EKa_' num2str(E_Ka) ...
'_epsilon_' num2str(epsilon) '_k_' num2str(k) '_n_' num2str(n) '.mat'];
else
filename = ['EbN0_SAMPR_KaPoissonKnown_EKa_' num2str(E_Ka) ...
'_epsilon_' num2str(epsilon) '_k_' num2str(k) '_n_' num2str(n) '.mat'];
end
save(filename, 'data', '-v7.3');
else
keyboard
end
end
%%
function PMF_Ksa = PMF_Ksa(p_Ka,E_Ka,L,Ksa)
% Compute the PMF of the number of active users per slot
PMF_Ksa = zeros(size(Ksa));
K_u = E_Ka;
while p_Ka(K_u+1) > eps
K_u = K_u + 1;
end
K_u = max(K_u,10000);
PMF_Ksa_scalar = @(T) sum(p_Ka(0:K_u).*binopdf(T,0:K_u,1/L));
for idxKsa = 1:length(Ksa)
PMF_Ksa(idxKsa) = PMF_Ksa_scalar(Ksa(idxKsa));
end
end
%%
function [rcu,P,P1,nSlot] = binary_search_P_SAMPR_KaKnown(f,x1,x2,TARGET,TOL,E_Ka)
% Search for the value of P such that the RCU bound of the PUPE of SAMPR
% given by min_{P1,L} f(P,P1,L) is between TARGET - TOL and TARGET
iter= 20; % maximum number of iterations
k1=0; % iteration index
[rcu_tmp,P1_tmp,nSlot_tmp] = golden_search_P1_SAMPR_KaKnown(f,1e-9,x2,x2/100,E_Ka);
[rcu_tmp1,P1_tmp1,nSlot_tmp1] = golden_search_P1_SAMPR_KaKnown(f,1e-9,x1,x1/100,E_Ka);
if x1 == x2
P = x1;
rcu = rcu_tmp;
P1 = P1_tmp;
nSlot = nSlot_tmp;
elseif rcu_tmp >= TARGET
warning('Impossible to achieve the target within the given range of parameter :( ');
P = x2;
rcu = rcu_tmp;
P1 = P1_tmp;
nSlot = nSlot_tmp;
elseif rcu_tmp1 <= TARGET
warning('All parameter values in the range can achieve the target :) ');
P = x1;
rcu = rcu_tmp1;
P1 = P1_tmp1;
nSlot = nSlot_tmp1;
else
[fx,P1_tmp,nSlot_tmp]=...
golden_search_P1_SAMPR_KaKnown(f,1e-9,(x1+x2)/2,(x1+x2)/200,E_Ka);
while (TARGET < fx || fx < TARGET - TOL) && (k1<iter || TARGET < fx)
if k1 > 0
[fx,P1_tmp,nSlot_tmp]=...
golden_search_P1_SAMPR_KaKnown(f,0,(x1+x2)/2,(x1+x2)/200,E_Ka);
end
if TARGET > fx
x2 = (x1+x2)/2; %set new end of interval
else
x1 = (x1+x2)/2; %replace as new start index
end
k1=k1+1;
end
rcu = fx;
P = x2;
P1 = P1_tmp;
nSlot = nSlot_tmp;
end
end
%%
function [rcu, P1,nSlot] = golden_search_P1_SAMPR_KaKnown(f, START_INT, END_INT, TOL, E_Ka)
% Minimize min_L f(P,P1,L) over P1 for given P
P = END_INT;
iter= 20; % maximum number of iterations
tau=double((sqrt(5)-1)/2); % golden proportion coefficient, around 0.618
k2=0; % number of iterations
x1=START_INT+(1-tau)*(END_INT-START_INT); % computing x values
x2=START_INT+tau*(END_INT-START_INT);
TOL_nSlot = 1;
nSlot_l = E_Ka;
nSlot_u = 6*E_Ka;
[f_x1,nSlot_1] = golden_search_nSlot_SAMPR_KaKnown(f,P,x1,nSlot_l,nSlot_u,TOL_nSlot); %f(P,x1); % % computing values in x points
[f_x2,nSlot_2] = golden_search_nSlot_SAMPR_KaKnown(f,P,x2,nSlot_l,nSlot_u,TOL_nSlot); %f(P,x2);
while ((abs(END_INT-START_INT)>TOL) && (k2<iter))
if f_x1 < f_x2 %function smaller in x1 than in x2
END_INT = x2; %make interval smaller
x2 = x1; %set new end of interval
x1 = START_INT+(1-tau)*(END_INT-START_INT); %find new beginning
f_x2 = f_x1;%already have value in x1
nSlot_2 = nSlot_1;
% [fMD_x1,fFA_x1] = f(P,x1,2*E_Ka);
[f_x1,nSlot_1] = ...
golden_search_nSlot_SAMPR_KaKnown(f, P,x1,nSlot_l,nSlot_u,TOL_nSlot); %f(P,x1); %%compute new value for new beginning
else
START_INT=x1; %function smaller in x2 than in x1 so set new start indx to x1
x1 = x2; %replace as new start index
x2=START_INT+tau*(END_INT-START_INT); %compute new end index
f_x1= f_x2;
nSlot_1 = nSlot_2;
% [fMD_x2,fFA_x2] = f(P,x2,2*E_Ka);
[f_x2,nSlot_2] = ...
golden_search_nSlot_SAMPR_KaKnown(f, P,x2,nSlot_l,nSlot_u,TOL_nSlot); %f(P,x2);
end
k2=k2+1;
end
if f_x1 < f_x2
P1=x1;
rcu = f_x1;
nSlot = nSlot_1;
else
P1=x2;
rcu = f_x2;
nSlot = nSlot_2;
end
end
%%
function [rcu, nSlot] = golden_search_nSlot_SAMPR_KaKnown(f, P,P1,START_INT,END_INT,TOL)
% Minimize f(P,P1,L) over L for given P and P1
iter= 20; % maximum number of iterations
tau=double((sqrt(5)-1)/2); % golden proportion coefficient, around 0.618
k2=0; % number of iterations
x1=max(round(START_INT+(1-tau)*(END_INT-START_INT)),0); % computing x values
x2=round(START_INT+tau*(END_INT-START_INT));
f_x1 = f(P,P1,x1); % computing values in x points
f_x2 = f(P,P1,x2);
while ((abs(END_INT-START_INT)>TOL) && (k2<iter))
if f_x1 < f_x2 %function smaller in x1 than in x2
END_INT = x2; %make interval smaller
x2 = x1; %set new end of interval
x1 = max(round(START_INT+(1-tau)*(END_INT-START_INT)),0); %find new beginning
f_x2 = f_x1;%already have value in x1
f_x1 = f(P,P1,x1); %compute new value for new beginning
else
START_INT=x1; %function smaller in x2 than in x1 so set new start indx to x1
x1 = x2; %replace as new start index
x2 = round(START_INT+tau*(END_INT-START_INT)); %compute new end index
f_x1= f_x2;
f_x2 = f(P,P1,x2); %
end
k2=k2+1;
end
if f_x1 < f_x2
nSlot=x1;
rcu = f_x1;
else
nSlot=x2;
rcu = f_x2;
end
end
EbN0_SAMPR_KaPoissonUnknown.m
function data = EbN0_SAMPR_KaPoissonUnknown(k, n, epsilon_MD, epsilon_FA, E_Ka, SlotIdxCoding)
% function data = EbN0_ALOHA_KaPoissonUnknown(k, n, epsilon, E_Ka, SlotIdxCoding)
% Find the minimal required EbN0 (in dB) such that the misdetection and
% false alarm probabilities achieved with slotted ALOHA (SA) with
% multi-packet reception (MPR) is below certain thresholds epsilon_MD and
% epsilon_FA, respectively, for a system with the number of active
% users following a Poisson distribution, and unknown. See Corollary 2 in
%
% [1] K.-H. Ngo, A. Lancho, G. Durisi and A. G. i. Amat, "Unsourced
% Multiple Access With Random User Activity," submitted to IEEE Trans. Inf.
% Theory, Jan. 2022.
%
% INPUTS
% k : number of bits per symbol
% n : framelength (number of complex DoFs)
% epsilon_MD : target misdetection probability
% epsilon_FA : target false alarm probability
% E_Ka : average number of active users
% SlotIdxCoding : 1 if slot-index coding is employed, 0 otherwise
%
% OUTPUT
% data : store the system parameters and the minimal required EbN0
addpath ../RCU_KaUnknown
tic
DEBUG = 0;
%% debugging mode
if DEBUG == 1
k = 128;
n = 19200;
epsilon_MD = .1;
epsilon_FA = .1;
E_Ka = 50;
SlotIdxCoding = 0;
end
%% Poissom PMF of Ka, can be modified to consider other distributions
p_Ka = @(K) poisspdf(K,E_Ka);
%% The range of power to search
EbN0db_lowest = -1;
EbN0db_highest = 16;
P_lowest = k*10^(EbN0db_lowest/10)/n;
P_highest = k*10^(EbN0db_highest/10)/n;
%% Function to compute the RCU bound
if SlotIdxCoding == 1
f_rcu = @(P,P1,L,DecRad) RCU_KaRandomUnknown(P*L,P1*L,DecRad,DecRad,k-floor(log2(L)),...
floor(n/L),E_Ka/L,@(Ksa) PMF_Ksa(p_Ka,E_Ka,L,Ksa));
else
f_rcu = @(P,P1,L,DecRad) RCU_KaRandomUnknown(P*L,P1*L,DecRad,DecRad,k,floor(n/L),...
E_Ka/L,@(Ksa) PMF_Ksa(p_Ka,E_Ka,L,Ksa));
end
%% Search for the minimal required EbN0
% The RCU is optimized over P1, L, and the decoding radius
[eps_RCU_MD, eps_RCU_FA, P_RCU,P1,L,DecRad] = binary_search_P_SAMPR(f_rcu, ...
P_lowest, P_highest,epsilon_MD,epsilon_FA,min(epsilon_MD,epsilon_FA)/100,E_Ka);
EbN0db_RCU = 10*log10(n*P_RCU/k);
%% Save the results
sim_time = toc;
data.EbN0db = EbN0db_RCU;
data.E_Ka = E_Ka;
data.p_Ka = 'Poisson';
data.eps_RCU_MD = eps_RCU_MD;
data.eps_RCU_FA = eps_RCU_FA;
data.epsilon_MD = epsilon_MD;
data.epsilon_FA = epsilon_FA;
data.DecRad = DecRad;
data.k = k;
data.n = n;
data.P1 = P1;
data.nSlot = L;
data.SlotIdxCoding = SlotIdxCoding;
data.sim_time = sim_time;
if DEBUG ~= 1
if SlotIdxCoding == 1
filename = ['EbN0_SAMPR_SlotIdxCoding_KaPoissonUnknown_EKa_' num2str(E_Ka)...
'_epsilonMD_' num2str(epsilon_MD) '_epsilonFA_' num2str(epsilon_FA)...
'_k_' num2str(k) '_n_' num2str(n) '.mat'];
else
filename = ['EbN0_SAMPR_KaPoissonUnknown_EKa_' num2str(E_Ka) ...
'_epsilonMD_' num2str(epsilon_MD) '_epsilonFA_' num2str(epsilon_FA) ...
'_k_' num2str(k) '_n_' num2str(n) '.mat'];
end
save(filename, 'data', '-v7.3');
else
keyboard
end
end
%%
function PMF_Ksa = PMF_Ksa(p_Ka,E_Ka,L,Ksa)
% Compute the PMF of the number of active users per slot
PMF_Ksa = zeros(size(Ksa));
K_u = E_Ka;
while p_Ka(K_u+1) > eps
K_u = K_u + 1;
end
K_u = max(K_u,10000);
PMF_Ksa_scalar = @(T) sum(p_Ka(0:K_u).*binopdf(T,0:K_u,1/L));
for idxKsa = 1:length(Ksa)
PMF_Ksa(idxKsa) = PMF_Ksa_scalar(Ksa(idxKsa));
end
end
%%
function [rcu_MD, rcu_FA, P,P1,nSlot,DecRad] = binary_search_P_SAMPR(f,x1,x2,TARGET_MD,TARGET_FA,TOL,E_Ka)
% Search for the value of P such that the RCU bound of the PUPE of SAMPR
% given by min_{P1,L,DecRad} weightedSum[f(P,P1,L,DecRad)] is between
% TARGET - TOL and TARGET
weight_MD = 1/(1 + TARGET_MD/TARGET_FA);
weight_FA = 1/(1 + TARGET_FA/TARGET_MD);
iter= 20; % maximum number of iterations
k1=0; % number of iterations
[rcu_MD_tmp,rcu_FA_tmp,P1_tmp,nSlot_tmp,DecRad_tmp] = ...
golden_search_P1_SAMPR(f,1e-9,x2,x2/100, weight_MD, weight_FA,E_Ka);
[rcu_MD_tmp1,rcu_FA_tmp1,P1_tmp1,nSlot_tmp1,DecRad_tmp1] = ...
golden_search_P1_SAMPR(f,1e-9,x1,x1/100, weight_MD, weight_FA,E_Ka);
if x1 == x2
P = x1;
rcu_MD = rcu_MD_tmp;
rcu_FA = rcu_FA_tmp;
P1 = P1_tmp;
nSlot = nSlot_tmp;
DecRad = DecRad_tmp;
elseif rcu_MD_tmp >= TARGET_MD || rcu_FA_tmp >= TARGET_FA
warning('Impossible to achieve the target within the given range of parameter :( ');
P = x2;
rcu_MD = rcu_MD_tmp;
rcu_FA = rcu_FA_tmp;
P1 = P1_tmp;
nSlot = nSlot_tmp;
DecRad = DecRad_tmp;
elseif rcu_MD_tmp1 <= TARGET_MD && rcu_FA_tmp1 <= TARGET_FA
warning('All parameter values in the range can achieve the target :) ');
P = x1;
rcu_MD = rcu_MD_tmp1;
rcu_FA = rcu_FA_tmp1;
P1 = P1_tmp1;
nSlot = nSlot_tmp1;
DecRad = DecRad_tmp1;
else
[fx_MD,fx_FA,P1_tmp,nSlot_tmp,DecRad_tmp]=...
golden_search_P1_SAMPR(f,1e-9,(x1+x2)/2,(x1+x2)/200, weight_MD, weight_FA,E_Ka);
while ~((TARGET_MD >= fx_MD && fx_MD >= TARGET_MD - TOL && TARGET_FA >= fx_FA) || ...
(TARGET_FA >= fx_FA && fx_FA >= TARGET_FA - TOL && TARGET_MD >= fx_MD) ...
|| (k1>iter))
if k1 > 0
[fx_MD,fx_FA,P1_tmp,nSlot_tmp,DecRad_tmp]=...
golden_search_P1_SAMPR(f,0,(x1+x2)/2,(x1+x2)/200, weight_MD, weight_FA,E_Ka);
end
if TARGET_MD > fx_MD && TARGET_FA > fx_FA
x2 = (x1+x2)/2; %set new end of interval
else
x1 = (x1+x2)/2; %replace as new start index
end
k1=k1+1;
end
rcu_MD = fx_MD;
rcu_FA = fx_FA;
P = x2;
P1 = P1_tmp;
nSlot = nSlot_tmp;
DecRad = DecRad_tmp;
end
end
%%
function [rcu_MD, rcu_FA, P1,nSlot,DecRad] = golden_search_P1_SAMPR(f, START_INT, END_INT, TOL, weight_MD, weight_FA,E_Ka)
% Optimize P1
P = END_INT;
iter= 20; % maximum number of iterations
tau=double((sqrt(5)-1)/2); % golden proportion coefficient, around 0.618
k2=0; % number of iterations
x1=START_INT+(1-tau)*(END_INT-START_INT); % computing x values
x2=START_INT+tau*(END_INT-START_INT);
% K_l = Ka;
% K_u = 1.5*E_Ka;
TOL_nSlot = 1;
nSlot_l = E_Ka;
nSlot_u = 6*E_Ka;
[fMD_x1,fFA_x1,nSlot_1,DecRad_1] = ...
golden_search_nSlot_SAMPR(f, P,x1,nSlot_l,nSlot_u,TOL_nSlot,weight_MD,weight_FA);
[fMD_x2,fFA_x2,nSlot_2,DecRad_2] = ...
golden_search_nSlot_SAMPR(f, P,x2,nSlot_l,nSlot_u,TOL_nSlot,weight_MD,weight_FA);
while ((abs(END_INT-START_INT)>TOL) && (k2<iter))
if fMD_x1*weight_MD + fFA_x1*weight_FA < fMD_x2*weight_MD + fFA_x2*weight_FA %function smaller in x1 than in x2
END_INT = x2; %make interval smaller
x2 = x1; %set new end of interval
x1 = START_INT+(1-tau)*(END_INT-START_INT); %find new beginning
fMD_x2 = fMD_x1;%already have value in x1
fFA_x2 = fFA_x1;
DecRad_2 = DecRad_1;
nSlot_2 = nSlot_1;
[fMD_x1,fFA_x1,nSlot_1,DecRad_1] = ...
golden_search_nSlot_SAMPR(f, P,x1,nSlot_l,nSlot_u,TOL_nSlot,weight_MD,weight_FA); %compute new value for new beginning
else
START_INT=x1; %function smaller in x2 than in x1 so set new start indx to x1
x1 = x2; %replace as new start index
x2=START_INT+tau*(END_INT-START_INT); %compute new end index
fMD_x1= fMD_x2;
fFA_x1 = fFA_x2;
nSlot_1 = nSlot_2;
DecRad_1 = DecRad_2;
[fMD_x2,fFA_x2,nSlot_2,DecRad_2] = ...
golden_search_nSlot_SAMPR(f, P,x2,nSlot_l,nSlot_u,TOL_nSlot,weight_MD,weight_FA);
end
k2=k2+1;
end
if fMD_x1*weight_MD + fFA_x1*weight_FA < fMD_x2*weight_MD + fFA_x2*weight_FA
P1=x1;
rcu_MD = fMD_x1;
rcu_FA = fFA_x1;
nSlot = nSlot_1;
DecRad = DecRad_1;
else
P1=x2;
rcu_MD = fMD_x2;
rcu_FA = fFA_x2;
nSlot = nSlot_2;
DecRad = DecRad_2;
end
end
%%
function [rcu_MD, rcu_FA, nSlot, DecRad] = golden_search_nSlot_SAMPR(f, P,P1,START_INT,END_INT,TOL,weight_MD,weight_FA)
% Optimize the number of slots
iter= 20; % maximum number of iterations
tau=double((sqrt(5)-1)/2); % golden proportion coefficient, around 0.618
k2=0; % number of iterations
x1=max(round(START_INT+(1-tau)*(END_INT-START_INT)),0); % computing x values
x2=round(START_INT+tau*(END_INT-START_INT));
[fMD_x1,fFA_x1,DecRad_x1] = linear_search_DecRad_SAMPR(f,P,P1,x1,weight_MD,weight_FA); % computing values in x points
[fMD_x2,fFA_x2,DecRad_x2] = linear_search_DecRad_SAMPR(f,P,P1,x2,weight_MD,weight_FA);
while ((abs(END_INT-START_INT)>TOL) && (k2<iter))
if fMD_x1*weight_MD + fFA_x1*weight_FA < fMD_x2*weight_MD + fFA_x2*weight_FA %function smaller in x1 than in x2
END_INT = x2; %make interval smaller
x2 = x1; %set new end of interval
x1 = max(round(START_INT+(1-tau)*(END_INT-START_INT)),0); %find new beginning
fMD_x2 = fMD_x1;%already have value in x1
fFA_x2 = fFA_x1;
DecRad_x2 = DecRad_x1;
if fixListSize
[fMD_x1,fFA_x1] = f(P,P1,x1); DecRad_x1 = 0;
else
[fMD_x1,fFA_x1,DecRad_x1] = linear_search_DecRad_SAMPR(f,P,P1,x1,weight_MD,weight_FA); %compute new value for new beginning
end
else
START_INT=x1; %function smaller in x2 than in x1 so set new start indx to x1
x1 = x2; %replace as new start index
x2 = round(START_INT+tau*(END_INT-START_INT)); %compute new end index
fMD_x1= fMD_x2;
fFA_x1 = fFA_x2;
DecRad_x2 = 0;
if fixListSize
[fMD_x2,fFA_x2] = f(P,P1,x2);
else
[fMD_x2,fFA_x2,DecRad_x2] = linear_search_DecRad_SAMPR(f,P,P1,x2,weight_MD,weight_FA);
end
end
k2=k2+1;
end
if fMD_x1*weight_MD + fFA_x1*weight_FA < fMD_x2*weight_MD + fFA_x2*weight_FA
nSlot=x1;
rcu_MD = fMD_x1;
rcu_FA = fFA_x1;
DecRad = DecRad_x1;
else
nSlot=x2;
rcu_MD = fMD_x2;
rcu_FA = fFA_x2;
DecRad = DecRad_x2;
end
end
%%
function [fMD,fFA,DecRad] = linear_search_DecRad_SAMPR(f,P,P1,x,weight_MD,weight_FA)
% Optimize the decoding radius
maxrad = 5; % search decoding radius from 0 to 4
fMD = zeros(maxrad,1);
fFA = zeros(maxrad,1);
for ii = 1:maxrad
[fMD_tmp, fFA_tmp] = f(P,P1,x,ii-1);
fMD(ii) = fMD_tmp;
fFA(ii) = fFA_tmp;
end
[~,idxMin] = min(fMD*weight_MD + fFA*weight_FA);
fMD = fMD(idxMin);
fFA = fFA(idxMin);
iiset = 0:maxrad-1;
DecRad = iiset(idxMin);
% DecRad = 0;
% [fMD, fFA] = f(P,P1,x,DecRad);
end
TIN
EbN0_TIN_KaPoissonUnknown.m
function data = EbN0_TIN_KaPoissonUnknown(k, n, epsilon_MD, epsilon_FA, E_Ka, normalApprox)
% function data = EbN0_TIN_KaPoissonUnknown(k, n, epsilon_MD, epsilon_FA, E_Ka, normalApprox)
% Find the minimal required EbN0 (in dB) such that the misdetection and
% false alarm probabilities achieved with treating-interference-as-noise
% (TIN) is below certain thresholds epsilon_MD and
% epsilon_FA, respectively, for a system with the number of active
% users following a Poisson distribution, and unknown. See Theorem 4 in
%
% [1] K.-H. Ngo, A. Lancho, G. Durisi and A. G. i. Amat, "Unsourced
% Multiple Access With Random User Activity," submitted to IEEE Trans. Inf.
% Theory, Jan. 2022.
%
% INPUTS
% k : number of bits per symbol
% n : framelength (number of complex DoFs)
% epsilon_MD : target misdetection probability
% epsilon_FA : target false alarm probability
% E_Ka : average number of active users
% normalApprox : 1 if normal approximation is used, 0 otherwise
%
% OUTPUT
% data : store the system parameters and the minimal required EbN0
tic
DEBUG = 0;
%% debugging mode
if DEBUG == 1
k = 128; % Number of bits
n = 19200; % Frame length
epsilon_MD = .1; % Per-user error probability
epsilon_FA = .1; % Per-user error probability
E_Ka = 100;
normalApprox = 1;
end
%% Poissom PMF of Ka, can be modified to consider other distributions
p_Ka = @(K) poisspdf(K,E_Ka);
%% The range of power to search
EbN0db_lowest = 0;
EbN0db_highest = 15;
P_lowest = k*10^(EbN0db_lowest/10)/n;
P_highest = k*10^(EbN0db_highest/10)/n;
%% Function to compute the RCUs-TIN bound
f_TIN = @(P,P1,s) RCUs_TIN(P,P1,s,k,n,E_Ka,p_Ka,normalApprox);
%% Search for the minimal required EbN0
[eps_TIN_MD, eps_TIN_FA, P,P1,sopt] = binary_search_TIN(f_TIN, P_lowest, P_highest,...
epsilon_MD,epsilon_FA,min(epsilon_MD,epsilon_FA)/100,normalApprox);
EbN0db = 10*log10(n*P/k);
%% Save the results
sim_time = toc;
data.EbN0db = EbN0db;
data.E_Ka = E_Ka;
data.p_Ka = 'Poisson';
data.eps_TIN_MD = eps_TIN_MD;
data.eps_TIN_FA = eps_TIN_FA;
data.epsilon_MD = epsilon_MD;
data.epsilon_FA = epsilon_FA;
data.k = k;
data.n = n;
data.P1 = P1;
data.s = sopt;
data.normalApprox = normalApprox;
data.sim_time = sim_time;
if DEBUG ~= 1
filename = ['EbN0_TIN_KaPoissonUnknown_EKa_' num2str(E_Ka) '_epsilonMD_' ...
num2str(epsilon_MD) '_epsilonFA_' num2str(epsilon_FA) ...
'_k_' num2str(k) '_n_' num2str(n) '.mat'];
save(filename, 'data', '-v7.3');
else
keyboard
end
end
%%
function Pe = eta_s(n,P,R,Ka_vec,s,N)
% Computation \eta_s in [1, th. 4], which is also a RCUs bound for the
% noncoherent Rayleigh block-fading channel.
%
% INPUTS:
% n = blocklength
% P = transmit power (no dB!)
% R = Values of Rate ln(M)/n at which the bound will be computed
% Ka = Number of active users
% s = parameter to optimize in the RCUs (s>0)
% N = Number of samples to generate the generalized info. density
%
% OUTPUT:
% Pe = Error probability obtained as a result of the computation of the bound
Pe = nan(size(Ka_vec));
parfor ii = 1:length(Ka_vec)
Ka = Ka_vec(ii);
snr = P/(1+(Ka-1)*P);
x = sqrt(0.5*snr)*(randn(n,N) + 1i*randn(n,N)); % signal tx by one of the users
z = sqrt(0.5)*(randn(n,N) + 1i*randn(n,N)); % Noise plus interference
y = x + z; % Received signal
i_s = -s*vecnorm(y-x).^2 + s*vecnorm(y).^2/(1+s*snr) + n*log(1+s*snr); % generalized information density
Pe(ii) = mean( exp( - max(0, i_s - log(exp(n*R)-1))));
end
end
%%
function [eps_MD,eps_FA] = RCUs_TIN(P,P1,s,k,n,E_Ka,p_Ka,normalApprox)
% codebook size
M = 2^k;
%% COMPUTATION OF \tilde{p}
K_l = floor(E_Ka); K_u = ceil(E_Ka);
while p_Ka(K_l-1) > 1e-9
K_l = K_l - 1;
end
K_l = max(K_l,0);
while p_Ka(K_u+1) > 1e-9
K_u = K_u + 1;
end
p01 = 0;
for Ka_tmp = K_l:K_u
p01_tmp = 1;
for ii = 1:Ka_tmp-1
p01_tmp = p01_tmp*(1-ii/M);
end
p01 = p01 + p_Ka(Ka_tmp)*p01_tmp;
end
ptilde = E_Ka*gammainc(n*P/P1,n,'upper') + 1 - p01 + 1 - sum(p_Ka(K_l:K_u));
%% Initialize eps_MD and eps_FA to be \tilde{p}
eps_MD = ptilde;
eps_FA = ptilde;
%% Compute the \xi(Ka,Ka') term for Ka and Ka' from K_l to K_u
if K_l == K_u
xi_Ka_KaEst = 1;
else
xi_Ka_KaEst = zeros(K_u-K_l+1,K_u-K_l+1);
for Ka = K_l:K_u
KaEst_thres = @(Ka,KaEst,P1) n.*(log(1+Ka.*P1) - log(1+KaEst.*P1))./((1+Ka.*P1)./(1+KaEst.*P1)-1);
P_Ka_KaEst_list_a = gammainc(KaEst_thres(Ka,K_l:Ka-1,P1),n,'lower');
P_Ka_KaEst_list_b = gammainc(KaEst_thres(Ka,Ka+1:K_u,P1),n,'upper');
P_Ka_Ka = 1 - max([P_Ka_KaEst_list_a P_Ka_KaEst_list_b]);
xi_Ka_KaEst(:,Ka-K_l+1) = [P_Ka_KaEst_list_a'; P_Ka_Ka; P_Ka_KaEst_list_b'];
end
end
%% Compute the \eta_s terms
N = 1e2;
R = k/n*log(2);
if normalApprox
C_TIN = @(Ka,P) log(1+P./(1+P.*(Ka-1)));
V_TIN = @(Ka,P) sqrt(2.*P.*(log(exp(1))).^2./(1+Ka.*P));
eta_s_vec = qfunc((n.*C_TIN((K_l:K_u)',P1) - log(M-1)) ./ (sqrt(n).*V_TIN((K_l:K_u)',P1)));
else
eta_s_vec = eta_s(n,P1,R,(K_l:K_u)',s,N);
end
%% The sum over Ka
parfor Ka = K_l:K_u
nAdditionalMDFA = min(min((K_l:K_u)',M-max(Ka,(K_l:K_u)')),Ka);
eta = min(eta_s_vec,xi_Ka_KaEst(Ka-K_l+1,:)');
% the sum over Ka' is implemented by the inner products
if Ka > 0
eps_MD = eps_MD + (feval(p_Ka,Ka)/Ka)*(nAdditionalMDFA'*eta + max(Ka-(K_l:K_u),0)*xi_Ka_KaEst(Ka-K_l+1,:)');
end
if K_l == 0
eps_FA = eps_FA + ...
feval(p_Ka,Ka)*sum((nAdditionalMDFA(2:end).*eta(2:end) + max((K_l+1:K_u)'-Ka,0).*xi_Ka_KaEst(Ka-K_l+1,2:end)')./(K_l+1:K_u)');
else
eps_FA = eps_FA + ...
feval(p_Ka,Ka)*sum((nAdditionalMDFA.*eta + max((K_l:K_u)'-Ka,0).*xi_Ka_KaEst(Ka-K_l+1,:)')./(K_l:K_u)');
end
% if eps_MD >= 1 && eps_FA >= 1
% break;
% end
end
eps_MD = min(eps_MD,1);
eps_FA = min(eps_FA,1);
end
%%
function [rcu_MD, rcu_FA, P,P1,s_opt] = binary_search_TIN(f,x1,x2,TARGET_MD,TARGET_FA,TOL,normalApprox)
% Search for P
weight_MD = 1/(1 + TARGET_MD/TARGET_FA);
weight_FA = 1/(1 + TARGET_FA/TARGET_MD);
iter= 20; % maximum number of iterations
k1=0; % number of iterations
[rcu_MD_tmp,rcu_FA_tmp,P1_tmp,s_tmp] = golden_search_P1_TIN(f,1e-9,x2,x2/200, weight_MD, weight_FA,normalApprox);
[rcu_MD_tmp1,rcu_FA_tmp1,P1_tmp1,s_tmp1] = golden_search_P1_TIN(f,1e-9,x1,x1/200, weight_MD, weight_FA,normalApprox);
if x1 == x2
P = x1;
rcu_MD = rcu_MD_tmp;
rcu_FA = rcu_FA_tmp;
P1 = P1_tmp;
s_opt = s_tmp;
elseif rcu_MD_tmp > TARGET_MD || rcu_FA_tmp > TARGET_FA
warning('Impossible to achieve the target within the given range of parameter :( ');
P = x2;
rcu_MD = rcu_MD_tmp;
rcu_FA = rcu_FA_tmp;
P1 = P1_tmp;
s_opt = s_tmp1;
elseif rcu_MD_tmp1 < TARGET_MD || rcu_FA_tmp1 < TARGET_FA
warning('All parameter values in the range can achieve the target :) ');
P = x1;
rcu_MD = rcu_MD_tmp1;
rcu_FA = rcu_FA_tmp1;
P1 = P1_tmp1;
s_opt = s_tmp1;
else
[fx_MD,fx_FA,P1_tmp,s_tmp]=golden_search_P1_TIN(f,1e-9,(x1+x2)/2,(x1+x2)/200, weight_MD, weight_FA,normalApprox); % computing values in x points
while ~((TARGET_MD >= fx_MD && fx_MD >= TARGET_MD - TOL && TARGET_FA >= fx_FA) || ...
(TARGET_FA >= fx_FA && fx_FA >= TARGET_FA - TOL && TARGET_MD >= fx_MD) ...
|| (k1>iter))
if k1 > 0
[fx_MD,fx_FA,P1_tmp,s_tmp]=golden_search_P1_TIN(f,0,(x1+x2)/2,(x1+x2)/200, weight_MD, weight_FA,normalApprox);
end
if TARGET_MD > fx_MD && TARGET_FA > fx_FA
x2 = (x1+x2)/2; %set new end of interval
else
x1 = (x1+x2)/2; %replace as new start index
end
k1=k1+1;
end
rcu_MD = fx_MD;
rcu_FA = fx_FA;
P = x2;
P1 = P1_tmp;
s_opt = s_tmp;
end
end
%%
function [rcu_MD, rcu_FA, P1,s_opt] = golden_search_P1_TIN(f, START_INT, END_INT, TOL, weight_MD, weight_FA,normalApprox)
% Optimize P1
P = END_INT;
iter= 20; % maximum number of iterations
tau=double((sqrt(5)-1)/2); % golden proportion coefficient, around 0.618
k2=0; % number of iterations
x1=START_INT+(1-tau)*(END_INT-START_INT); % computing x values
x2=START_INT+tau*(END_INT-START_INT);
s_start = 0;
s_end = 2;
if normalApprox
[fMD_x1,fFA_x1] = f(P,x1,1); s_opt1 = 1;
[fMD_x2,fFA_x2] = f(P,x2,1); s_opt2 = 1;
else
[fMD_x1,fFA_x1,s_opt1] = golden_search_s_TIN(f, P, x1, s_start, s_end, TOL, weight_MD, weight_FA);
[fMD_x2,fFA_x2,s_opt2] = golden_search_s_TIN(f, P, x2, s_start, s_end, TOL, weight_MD, weight_FA);
end
while ((abs(END_INT-START_INT)>TOL) && (k2<iter))
if fMD_x1*weight_MD + fFA_x1*weight_FA < fMD_x2*weight_MD + fFA_x2*weight_FA %function smaller in x1 than in x2
END_INT = x2; %make interval smaller
x2 = x1; %set new end of interval
x1 = START_INT+(1-tau)*(END_INT-START_INT); %find new beginning
fMD_x2 = fMD_x1; %already have value in x1
fFA_x2 = fFA_x1;
s_opt2 = s_opt1;
if normalApprox
[fMD_x1,fFA_x1] = f(P,x1,1); s_opt1 = 1;
else
[fMD_x1,fFA_x1,s_opt1] = golden_search_s_TIN(f, P, x1, s_start, s_end, TOL, weight_MD, weight_FA); %%compute new value for new beginning
end
else
START_INT=x1; %function smaller in x2 than in x1 so set new start indx to x1
x1 = x2; %replace as new start index
x2=START_INT+tau*(END_INT-START_INT); %compute new end index
fMD_x1= fMD_x2;
fFA_x1 = fFA_x2;
s_opt1 = s_opt2;
if normalApprox
[fMD_x2,fFA_x2] = f(P,x2,1); s_opt2 = 1;
else
[fMD_x2,fFA_x2,s_opt2] = golden_search_s_TIN(f, P, x2, s_start, s_end, TOL, weight_MD, weight_FA);
end
end
k2=k2+1;
end
if fMD_x1*weight_MD + fFA_x1*weight_FA < fMD_x2*weight_MD + fFA_x2*weight_FA
P1=x1;
s_opt = s_opt1;
rcu_MD = fMD_x1;
rcu_FA = fFA_x1;
else
P1=x2;
s_opt = s_opt2;
rcu_MD = fMD_x2;
rcu_FA = fFA_x2;
end
end
%%
function [rcu_MD, rcu_FA, s_opt] = golden_search_s_TIN(f, P, P1, START_INT, END_INT, TOL, weight_MD, weight_FA)
% Optimize s
iter= 20; % maximum number of iterations
tau=double((sqrt(5)-1)/2); % golden proportion coefficient, around 0.618
k2=0; % number of iterations
x1=START_INT+(1-tau)*(END_INT-START_INT); % computing x values
x2=START_INT+tau*(END_INT-START_INT);
[fMD_x1,fFA_x1] = f(P,P1,x1);
[fMD_x2,fFA_x2] = f(P,P1,x2);
while ((abs(END_INT-START_INT)>TOL) && (k2<iter))
if fMD_x1*weight_MD + fFA_x1*weight_FA < fMD_x2*weight_MD + fFA_x2*weight_FA %function smaller in x1 than in x2
END_INT = x2; %make interval smaller
x2 = x1; %set new end of interval
x1 = START_INT+(1-tau)*(END_INT-START_INT); %find new beginning
fMD_x2 = fMD_x1;%already have value in x1
fFA_x2 = fFA_x1;
[fMD_x1,fFA_x1] = f(P,P1,x1); %%compute new value for new beginning
else
START_INT=x1; %function smaller in x2 than in x1 so set new start indx to x1
x1 = x2; %replace as new start index
x2=START_INT+tau*(END_INT-START_INT); %compute new end index
fMD_x1= fMD_x2;
fFA_x1 = fFA_x2;
[fMD_x2,fFA_x2] = f(P,P1,x2);
end
k2=k2+1;
end
if fMD_x1*weight_MD + fFA_x1*weight_FA < fMD_x2*weight_MD + fFA_x2*weight_FA
s_opt=x1;
rcu_MD = fMD_x1;
rcu_FA = fFA_x1;
else
s_opt=x2;
rcu_MD = fMD_x2;
rcu_FA = fFA_x2;
end
end