仿真内容

文件中包含了一个关于无源多用户接入（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) 中展示。

如果你使用了这些代码，请引用上述的论文。

代码介绍部分

1. /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

 辅助函数。

2. /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

 辅助函数。

3. /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]。

4. /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 = k10^(EbN0db_lowest/10)/n;

P_highest = k10^(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 = 10log10(nP_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./nlog(M)-1./(nt).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 + 4P1.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(-nEt_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_Kagammainc(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) fMDweight_MD + fFAweight_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 &gt;= fx_MD &amp;&amp; fx_MD &gt;= TARGET_MD - TOL &amp;&amp; TARGET_FA &gt;= fx_FA) || ...
        (TARGET_FA &gt;= fx_FA &amp;&amp; fx_FA &gt;= TARGET_FA - TOL &amp;&amp; TARGET_MD &gt;= fx_MD) ...
        || (k1&gt;iter)) 
    if k1 &gt; 0
        [fx_MD,fx_FA,P1_tmp] = search_P1((x1+x2)/2); 
    end
    if TARGET_MD &gt; fx_MD &amp;&amp; TARGET_FA &gt; 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 = k10^(EbN0db_lowest/10)/n;

P_highest = k10^(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 = 10log10(nP_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(.0001min(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(.0001min(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_Kagammainc(nP/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 = k10^(EbN0db_lowest/10)/n;

P_highest = k10^(EbN0db_highest/10)/n;
%% Function to compute the RCU bound

if SlotIdxCoding == 1

f_rcu = @(P,P1,L) RCU_KaRandomKnown(PL,P1L,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(PL,P1L,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 = 10log10(nP_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,2E_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,2E_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 = k10^(EbN0db_lowest/10)/n;

P_highest = k10^(EbN0db_highest/10)/n;
%% Function to compute the RCU bound

if SlotIdxCoding == 1

f_rcu = @(P,P1,L,DecRad) RCU_KaRandomUnknown(PL,P1L,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(PL,P1L,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 = 10log10(nP_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_x1weight_MD + fFA_x1weight_FA < fMD_x2weight_MD + fFA_x2weight_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_x1weight_MD + fFA_x1weight_FA < fMD_x2weight_MD + fFA_x2weight_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_x1weight_MD + fFA_x1weight_FA < fMD_x2weight_MD + fFA_x2weight_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_x1weight_MD + fFA_x1weight_FA < fMD_x2weight_MD + fFA_x2weight_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(fMDweight_MD + fFAweight_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 = k10^(EbN0db_lowest/10)/n;

P_highest = k10^(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 = 10log10(nP/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.5snr)(randn(n,N) + 1irandn(n,N)); % signal tx by one of the users

z   = sqrt(0.5)(randn(n,N) + 1irandn(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_Kagammainc(nP/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_x1weight_MD + fFA_x1weight_FA < fMD_x2weight_MD + fFA_x2weight_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_x1weight_MD + fFA_x1weight_FA < fMD_x2weight_MD + fFA_x2weight_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_x1weight_MD + fFA_x1weight_FA < fMD_x2weight_MD + fFA_x2weight_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_x1weight_MD + fFA_x1weight_FA < fMD_x2weight_MD + fFA_x2weight_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

本栏目推荐文章
• ACCESS 快速构建修改数据的窗体
• 解决前后端的跨域问题：Access to XMLHttpRequest at '**' from origin '**' has been blocked by CORS policy: No 'Access-Control-Allow-Origin' header is present on the requested resource.
• android编译kanzi 问题 (1) Caused by: java.io.FileNotFoundException: Error: Could not find or access Kanzi's Android libraries directory: C:\XXXXXXX
• [LeetCode] 1363. Largest Multiple of Three 形成三的最大倍数
• No 'Access-Control-Allow-Origin' header is present on the requested resource'， 跨域访问的解决方法
• Android——2个activity之间的数据传递返回
• 【常用模块及方法】random_json_os_hashlib_time
• 【常用模块及方法】random_json_os_hashlib_time
• Jax框架在不同版本间的随机数生成并不一致 —— jax.random
• Android 优雅的Activity回调代码封装

Unsourced Multiple Activity Access Randomunsourced multiple activity access regression crc-aided unsourced access access-control-allow-origin 39 contains multiple unsourced multiple activity multiple-select activity7 activity_main articles learning multiple insomnia