Unsourced Multiple Access With Random User Activity
一、仿真内容
随记中包含了一个关于无源多用户接入（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) 中展示。 *** `如果你使用了这些代码，请引用上述的论文。`
m文件目录及介绍
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 的研究组提供，因此未包含在此随记中。
**注意**:每段代码中所用到的相关概念和MATLAB函数都附在了代码块的后面，方便查找。
二、RCU_KaKnown
2.1 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
%这个函数用于在给定目标值 TARGET 和容忍误差 TOL 的情况下，通过二分搜索算法在[x1 ，x2]区间中找到函数 f 的一个输入 P，使得 min_{P1} f(P, P1) 的值在区间 (TARGET - TOL, TARGET) 内。
iter = 20; % 最大迭代次数
k1 = 0; % 迭代次数计数
% 使用黄金分割算法计算 min_{P1} f(P,P1)
fx = golden_search(f,1e-9,(x1+x2)/2,(x1+x2)/200);
%调用了另一个函数 golden_search(见下)
while (TARGET < fx || fx<TARGET - TOL) && (k1<iter || TARGET < fx)
%fx的值在TARGET-TOL到TARGET或者达到最大迭代次数时结束循环，否则一直循环
if k1 > 0
fx = golden_search(f,0,(x1+x2)/2,(x1+x2)/200);
%第一次迭代后，重新使用黄金分割算法计算 min_{P1} f(P, P1)
end
if TARGET > fx %TARGET相当于f1，fx相当于f2
x2 = (x1+x2)/2; %迭代后若f2小于f1，证明最优值在前半段
else
x1 = (x1+x2)/2; %迭代后若f1小于f2，证明最优值在后半段
end
k1=k1+1; %迭代次数+1
end
epsilon = fx;%最终返回fx的值为epsilon
P = x1;%返回最优解，P就是使得 min_{P1} f(P, P1) 落在 (TARGET - TOL, TARGET) 区间内的结果。
end
2.2 golden_search.m
function [epsilon, P1] = golden_search(f, START_INT, END_INT, TOL)
% minimize f(P1) over P1. epsilon = min_{P1} f(P1)
%在区间[START_INT, END_INT]中寻找一个局部最优解P1，最小化f(P1)，搜索容忍误差为TOL
P = END_INT;
iter = 20; %最大迭代次数
tau = double((sqrt(5)-1)/2); % 黄金比例系数, 大约 0.618
k2 = 0; % 迭代次数计数
x1 = START_INT+(1-tau)*(END_INT-START_INT);%x1=a+0.382*（b-a）
x2 = START_INT+tau*(END_INT-START_INT);%x2=a+0.618*（b-a）
f_x1 = f(P,x1); %计算x1 x2起始处的函数值
f_x2 = f(P,x2);
while ((abs(END_INT-START_INT)>TOL) && (k2<iter))
%在b-a<TOL时或者达到最大迭代次数时才结束循环，否则一直循环
if f_x1 < f_x2 %f(x1)<f(x2) 最优解在(a,x2)之间
END_INT = x2; % x2处作为新区间终点
x2 = x1;
x1 = START_INT+(1-tau)*(END_INT-START_INT); %新的x1=a+0.382（x2-a）
f_x2 = f_x1;
f_x1 = f(P,x1);
else %f(x1)>=f(x2) 最优解在(x1,b)之间
START_INT=x1; % 从x1开始
x1 = x2;
x2=START_INT+tau*(END_INT-START_INT); %新的x2=x1+0.618（b-x1）
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
2.3 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
%该函数用于找到使 PUPE 低于给定阈值 epsilon 的最小所需 EbN0 值(每比特能量/噪声功率)
%针对的场景是Ka服从泊松分布且已知的情况
%输入参数：接收符号比特数k、帧长n、目标PUPE epsilon、平均活跃用户数 E_Ka
tic %tic函数用于记录函数执行的时间
DEBUG = 1;%DEBUG用于控制是否使用调试模式
%% 调试模式下的默认参数
if DEBUG == 1
k = 100;
n = 15000;
epsilon = 0.1;
E_Ka = 2;
end
%% 定义一个匿名函数 p_Ka ,表示活跃用户服从泊松分布的概率质量函数(PMF)
p_Ka = @(K) poisspdf(K,E_Ka); %使用MATLAB中的poisspdf函数计算泊松分布的概率，该函数见下,K和E_Ka时输入参数，E_Ka是给定时间间隔内事件发生的平均次数，函数返回K处计算的泊松分布pdf值
%% 设置搜索最小和最大EbN0值得初始范围
EbN0db_lowest = -0.5;%单位为dB
EbN0db_highest = 3;
P_lowest = k*10^(EbN0db_lowest/10)/n;%10^(EbN0db_lowest/10)计算的是每比特能量/噪声功率，dB转为比值，一个符号有kbit，帧长为n，所以要想真正计算每个bit能量与噪声比值需乘以k/n
P_highest = k*10^(EbN0db_highest/10)/n;
%% 定义匿名函数 f_rcu 来计算系统的 RCU（Random Coding Unbounded）界限
f_rcu = @(P,P1) RCU_KaRandomKnown(P,P1,k,n,E_Ka,p_Ka);%用到了RCU_KaRandomKnown函数(见下一个.m代码)
%% 使用二分搜索算法找到满足 PUPE小于epsilon 的最小EbN0值
[eps_RCU, P_RCU] = binary_search(f_rcu, P_lowest, P_highest, epsilon,epsilon/100);
EbN0db_RCU = 10*log10(n*P_RCU/(k));%乘以n/k再转为dB即为所求的最小EbN0的值
%[epsilon, P] = binary_search(f,x1,x2,TARGET,TOL)，其中TOL为epsilon/100
%% 计算算法执行时间以及结果参数，并存储到data结构体中
sim_time = toc;%tic是开始时间，记为0，toc是结束时间，所以toc就是运行时间
data.EbN0db = EbN0db_RCU;%matlab中不需要定义结构体直接使用就可以
data.E_Ka = E_Ka;
data.eps_est= eps_RCU;%最优值，就是把最后的f_rcu的值赋给了eps_RCU
data.epsilon= epsilon;%目标值
if DEBUG ~= 1%不是调试模式，将结果保存在文件中
filename = ['EbN0_KaPoissonKnown_EKa_' num2str(E_Ka) '_epsilon_'
num2str(epsilon) '_k_' num2str(k) '_n_' num2str(n) '.mat'];
%文件名基于输入参数E_Ka，epsilon，k，n，以便于标识不同的实验条件
save(filename, 'data', '-v7.3');%将data结构体保存到指定文件中
%'-v7.3' 表示以MATLAB 7.3版本的格式保存文件，这是一种支持大型数据的格式。
else
keyboard %如果是处于调试模式，程序暂停并等待键盘输入
end
end
poisspdf函数： 泊松概率密度
2.4 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])
%这段代码是用于计算Ka服从泊松分布且未知的情况下，PUPE和RCU界限
%输入参数为k，n，E_Ka和每比特能量EbN0db，输出存入结构体data
tic%计时
DEBUG = 0;%调试开关
%% 调试模式下默认参数
if DEBUG == 1
k = 128;
n = 19200;
E_Ka = 50;
EbN0db = 0;
end
%% 计算不同EbN0值对应的符号功率预算P_list
%注：功率预算是一种用于评估和分配通信系统中不同部分之间功率需求的概念。
P_list = k.*10.^(EbN0db./10)./n;%P_list是一个向量，表示真正的每个bit能量与噪声比值，其中每个元素对应一个不同的EbN0db
%% 初始化两个矩阵(设为全0)，用于存储计算结果
epsilon = zeros(length(EbN0db),length(E_Ka));
P1 = zeros(length(EbN0db),length(E_Ka));
%% 计算RCU界限
for idxEKa = 1:length(E_Ka)%嵌套循环
E_Ka = E_Ka(idxEKa);%更新活跃用户平均数E_Ka，将其设为当前迭代中的值，以便内层循环使用
p_Ka = @(K) poisspdf(K,E_Ka);%在E_Ka条件下，随机变量K的概率
for idxEbN0 = 1:length(EbN0db)%遍历EbN0db中的每一个元素
P = P_list(idxEbN0);%从P_list中选择当前循环迭代中的值，作为符号功率预算
f_rcu = @(P,P1) RCU_KaRandomKnown(P,P1,k,n,E_Ka,p_Ka);%定义匿名函数来计算RCU，P和P1分别代表符号功率预算和实际使用的功率，还使用了外层的E_Ka和p_Ka，返回RCU到f_rcu
[Pe_tmp,P1_tmp] = golden_search(f_rcu,0,P,P/100); % 使用黄金分割算法，返回最优解P1_tmp和最优值Pe_tmp
%将计算的最优值和最优解存入epsilon和P1矩阵中
epsilon(idxEbN0,idxEKa) = Pe_tmp;
P1(idxEbN0,idxEKa) = P1_tmp;
end
end
%% 存储结果到data结构体
sim_time = toc;%运行时间
data.EbN0db = EbN0db;
data.E_Ka = E_Ka;
data.p_Ka = 'Poisson';%表示使用的概率分布为泊松分布
data.epsilon= epsilon;%最优值矩阵（RCU界限矩阵）
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'); %见上一个.m
else
keyboard
end
end
2.5 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.
% 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;
%这段代码是用于计算Ka固定不变且已知的情况下的RCU界限
%输入参数：k，n，Ka，EbN0db
%输出参数：界限值epsilon、符号功率预算P、最优功率P1_opt
%% 符号功率预算P和码本大小M
P = k*10^(EbN0db/10)/n;
M = 2^k;
%% 初始化
% for the sum over t in Eq. (3)
t_vec = 1:Ka;%创建向量t_vec包含1到Ka的整数，用于后续计算
% for the optimization over rho and rho1 in E(t)
rho_vec = linspace(0,1,100); % 生成一个包含100个均匀分布在[0,1]的数值向量
rho1_vec= linspace(0,1,100);
rho = repmat(rho_vec,length(rho1_vec),1,length(t_vec));
%使用 repmat 函数把 rho_vec 复制成一个矩阵 rho ，length(rho1_vec) 表示在第一个维度上重复 rho_vec 的次数，1 表示在第二个维度上不进行重复，length(t_vec) 表示在第三个维度上重复 rho_vec 的次数。这样生成的 rho 矩阵是一个三维矩阵，其维度分别对应于 length(rho1_vec)、1、length(t_vec)
rho1 = permute(repmat(rho1_vec,length(rho_vec),1,length(t_vec)),[2,1,3]);
%首先使用 repmat函数 把 rho1_vec 复制成一个矩阵，length(rho1_vec) 表示在第一个维度上重复 rho_vec 的次数，1 表示在第二个维度上不进行重复，length(t_vec) 表示在第三个维度上重复 rho_vec 的次数。这样生成的 rho 矩阵是一个三维矩阵，其维度分别对应于 length(rho1_vec)、1、length(t_vec)。
% permute(..., [2, 1, 3]) 对生成的矩阵进行维度交换，将第一个和第二个维度交换。最终生成的 rho1 矩阵也是一个三维矩阵，其维度分别对应于 1、length(rho1_vec)、length(t_vec)。
t=permute(repmat(t_vec,length(rho1_vec),1,length(rho_vec)),[1,3,2]);
%同样，是在三个维度上重复 t_vec 向量，生成一个三维矩阵，其中维度分别对应于 length(rho1_vec)、1、length(rho_vec)。
% permute(..., [1, 3, 2]) 对生成的矩阵进行维度交换，将第二个和第三个维度交换。最终生成的 t 矩阵也是一个三维矩阵，其维度分别对应于 length(rho1_vec)、length(rho_vec)、1。
% Some functions for the computation of pt as defined in [1, Th. 1]
%定义了一系列函数用于计算RCU界限中的各种参数
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));%计算pt
% number of samples for empirical evaluation of the CDF of It in qt
%定义了用于计算 qt 的采样数 Niq ，该变量表示计算累积分布函数（CDF）时的采样数量
Niq = 1000;
% for the optimization over P' (here denoted by P1)
P1_vec = linspace(1e-9,P,20);%20各元素，在[1e-9，p]之间均匀分布
epsilon = zeros(size(P1_vec));% epsilon 初始化为一个全零向量，用于存储每个 P1 对应的 PUPE。
%% 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
%计算PUPE公式中的一项p0，它表示码本中两个符号的碰撞概率
p0 = nchoosek(Ka,2)/M + Ka*gammainc(n*P/P1,n,'upper');
%调用了nchoosek函数，表示从Ka元素中选择2个元素的组合数
%nchoosek(Ka,2)/M 表示在码本中两个符号发生碰撞的概率
%Ka*gammainc(n*P/P1,n,'upper')表示码本中每个符号发生错误的概率
%gammainc 是不完全伽马函数，'upper'表示计算的是上不完全伽马函数
%最终，p0 表示在码本中的两个符号中发生碰撞的概率与每个符号发生错误的概率之和。这是 PUPE 公式中的一个关键组成部分，用于描述整个码本的错误概率。
% Computation of pt
%计算PUPE公式中的另一项pt，它表示码本中每个符号的误码概率
pt = pt_f(P1,t,rho,rho1,n,M,Ka);
% Computation of qt (for t = 1 only, as in [1])
%计算PUPE公式中的另一项qt，它表示码本中第一个符号的误码概率，其中t=1
It = zeros(1,1000);
for II = 1:Niq %生成复高斯随机变量和随机码本
Zi = sqrt(.5)*(randn(1,n) + 1i*randn(1,n)); %生成大小为1*N的复高斯随机向量
codebook = sqrt(.5*P1)*(randn(Ka,n) + 1i*randn(Ka,n));%生成大小为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));%计算 t=1 时的统计量，其中使用了复数运算和矩阵重复。
It(II) = min(it); %记录每次迭代中的最小统计量
end
[prob,gam] = ecdf(It); %计算统计量 It 的经验累积分布函数（ECDF），结果存储在 prob 和 gam 中，调用的ecdf函数见代码块下
qt = min(prob + exp(n*(R1_f(n,M,1)+R2_f(n,Ka,1))-gam));
%计算 PUPE 公式中的 qt，其中 t=1。这涉及对 ECDF 进行一些调整和指数运算
% Computation of RHS of [1, Eq. (3)] for a given P' < P
% 在给定的P' < P 的情况，计算PUPE公式[1, Eq. (3)] 中的右侧项
epsilon(pp) = t_vec(1)/Ka*min(pt(1),qt) + t_vec(2:end)/Ka*pt(2:end) + p0;
%t_vec(1)/Ka*min(pt(1),qt): 计算 t=1 时的部分，其中 min(pt(1),qt) 表示取 pt(1) 和 qt 中的较小值。
%t_vec(2:end)/Ka*pt(2:end): 计算 t 大于等于 2 时的部分，其中 pt(2:end) 表示 pt 中从第二个元素开始的所有元素。
end
% Take the minimum over P'
%在所有计算的PUPE值中，找到最小值，并记录最优的P'
[epsilon,idx_P1opt] = min([epsilon 1]); %PUPE一定小于1，最小值存入epsilon，最小值的下标计入idx_P1opt
P1_opt = P1_vec(idx_P1opt);%通过下标找到epsilon对应的P'(即P1)，存入P1_opt
cdf(累积分布函数)：是一个概率分布的函数，表示随机变量小于或等于给定值的概率。对于一个随机变量 X 和一个实数 x，CDF(X) 给出了 X 小于或等于 x 的概率。
repmat函数：重复数组副本
permute 函数：置换数组维度
gammaln函数：gamma函数的对数
gammainc函数：不完全gamma函数
具体什么是上/下不完全伽马函数可以去百度。
nchoosek函数：二项式系数或所有组合
也就是我们学的 Ck n
ecdf函数(经验累积分布函数)：
[f，x]=ecdf(y)使用向量y中的数据返回经验累积分布函数(cdf)f，在x的各个点上计算。
2.6 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. Ka是随机未知的
% 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，即Ka的PMF
%
% 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
%用于对 It 的累积分布函数（CDF）进行经验估计的样本数量。
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 和 K_u ，使Ka 不在区间[K_l:K_u] 的概率很小。
K_l = floor(E_Ka); K_u = ceil(E_Ka);
%K_l 初始化为 E_Ka 向下取整，K_u 初始化为 E_Ka 向上取整。
while p_Ka(K_l-1) > 1e-9
%缩小 K_l 直到 p_Ka(K_l-1) 小于等于 1e-9，即落在该区间外的概率足够小。
K_l = K_l - 1;
end
K_l = max(K_l,0); %保证K_l>=0
while p_Ka(K_u+1) > 1e-9
%同样，缩小 K_u 直到 p_Ka(K_u+1) 小于等于 1e-9。
K_u = K_u + 1;
end
%计算 p0，其中包含了一些额外的项，考虑了 Ka 的随机性。
p01 = 0;%用于存储循环中计算得到的值
for Ka_tmp = K_l:K_u %嵌套循环，外层循环从 K_l 遍历到 K_u
p01_tmp = 1;
for ii = 1:Ka_tmp-1 %内层循环，计算每个 Ka_tmp 对应的一项 p01_tmp
p01_tmp = p01_tmp*(1-ii/M);
end
p01 = p01 + p_Ka(Ka_tmp)*p01_tmp;%将每个p01_tmp与相应的p_Ka(Ka_tmp)相乘，累加到p01
% 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');
% sum(p_Ka(K_l:K_u)) 表示在区间 [K_l, K_u] 内的 p_Ka 值(落入该区间的总概率)的和
%E_Ka*gammainc(n*P/P1,n,'upper')表示额外项
%% Initialize epsilon to p0
epsilon = p0;
%% Compute the RCU bound, averaged over the distribution of Ka
%对于每个Ka，并行计算 pt 和 qt
parfor Ka = max(K_l,1):K_u %parfor是并行循环，意思是每个 Ka 的迭代是独立的，因此可以并行执行，提高整体计算速度
%% Computation of pt
% Initialize for the sum over t and the optimization over rho and rho1
%在不同的 rho（符号能量分配系数）和 rho1（干扰能量分配系数）的情况下，计算了 pt(概率密度函数) 和 qt(分布函数)，并取二者第一个元素(t=1)的最小值来计算了 PUPE 边界 epsilon。
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]
%定义一系列匿名函数来计算pt
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 = pt_f(P1,t,rho,rho1,n,M,Ka);
%% Computation of qt (for t = 1 only, as in [1])，仅计算t=1时的qt
qt = 1;
% compute qt for Ka < 50 only due to complexity issue
if Ka < 50 %仅在Ka<50的条件下计算qt
It = zeros(1,Niq);
for II = 1:Niq %通过 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);%通过循环生成It
end
[prob,gam] = ecdf(It);%prob 是 It 的经验累积分布函数的概率
qt = min(prob + exp(n*(R1_f(n,M,1)+R2_f(n,Ka,1))-gam));%计算qt
end
%% Take the min of pt and qt
pt(1) = min(pt(1),qt);取pt(1)%与qt的最小值作为最终计算PUPE边界的参数
%% Compute epsilon, the RHS of [1, Eq. (3)]
epsilon = epsilon + (t_vec/Ka*pt)*feval(p_Ka,Ka);
end
end
三、RCU_KaUnknown
3.1 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.
%二分搜索算法，寻找使得 eps_MD 在 [TARGET_MD - TOL, TARGET_MD] 范围内，
% 同时 eps_FA 在 [TARGET_FA - TOL, TARGET_FA] 范围内的 P 的值。
% 其中 eps_MD 和 eps_FA 是通过调用函数 f(P, P1) 计算得到的，
% f 是在文献 [1] 中给出的 Theorem 1 中的 RCU 界限。
iter= 20; % 最大迭代次数
k1=0; % 目前迭代次数
%% 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; % P1 的固定比例，即 P1 = fracP1 * P，提高搜索效率
if fixP1 == 0 % 根据是否固定 P1 选择不同的搜索策略，等于0是不固定
if TARGET_MD == TARGET_FA % 如果目标 MD 和 FA 相等，则使用黄金分割法搜索 P1
search_P1 = @(x) golden_search_P1_MDFA(f,1e-9,x,x/100,'max',[],[]);
else % 如果目标 MD 和 FA 不相等，则计算基于目标概率的权重
weight_MD = 1/(1 + TARGET_MD/TARGET_FA);
weight_FA = 1/(1 + TARGET_FA/TARGET_MD);
% 使用带权重的黄金分割法搜索 P1
search_P1 = @(x) golden_search_P1_MDFA(f,1e-9,x,x/100,'weighted',weight_MD, weight_FA);
end
else %固定P1，则直接计算 MD 和 FA
search_P1 = @(x) compute_MDFA_fixedP1(f, x, fracP1); %函数定义见下
end
%% 寻找P
[eps_MD_tmp,eps_FA_tmp,P1_tmp] = search_P1(x2);%使用 search_P1 函数，传入 x2 作为起始 P 值，得到 eps_MD（MD 的界限）、eps_FA（FA 的界限）和 P1（在给定 P 下的最优 P1 值）。
[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; % 如果在给定功率范围内无法达到目标，则选择 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; % 如果所有功率值都能够达到目标，则选择 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))%终止条件是达到了指定的迭代次数或者满足了目标条件，即 eps_MD 在目标范围内，并且 eps_FA 在目标范围内
if k1 > 0
[fx_MD,fx_FA,P1_tmp] = search_P1((x1+x2)/2); %再次在 P 的当前区间中间位置进行一次搜索，获取新的 eps_MD、eps_FA 和 P1
end
if TARGET_MD > fx_MD && TARGET_FA > fx_FA
x2 = (x1+x2)/2; % 设置新的区间结束值
else
x1 = (x1+x2)/2; % 设置新的区间起始值
end
k1 = k1+1;
end
%返回最终结果
eps_MD = fx_MD;
eps_FA = fx_FA;
P = x2;
P1 = P1_tmp;
end
end
%定义函数compute_MDFA_fixedP1，计算固定P1时的 eps_MD, eps_FA, P1
function [eps_MD, eps_FA, P1] = compute_MDFA_fixedP1(f, P, fracP1)
P1 = P*fracP1;
[eps_MD, eps_FA] = f(P,P1);
end
3.2 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.
% 函数的目标是在区间 [0, P] 上找到一个最小值P1，这个最小值是由另一个函数 f(P,P1) 返回的 eps_MD 和 eps_FA 的权重和。这个权重和由参数 type 决定
if strcmp(type,'max') %不同类型g取值不同
g = @(fMD,fFA) max(fMD,fFA);
elseif strcmp(type,'weighted')
g = @(fMD,fFA) fMD*weight_MD + fFA*weight_FA;
end
P = END_INT; % 初始化功率预算P为结束边界END_INT
tau = double((sqrt(5)-1)/2); % 黄金分割法的比例参数 tau，大约0.618
iter = 20;
k2 = 0;
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) %g(x1)<g(x2)
END_INT = x2; %改变区间长度
x2 = x1;
x1 = START_INT+(1-tau)*(END_INT-START_INT);
fMD_x2 = fMD_x1;
fFA_x2 = fFA_x1;%上一轮循环中已经计算过了fMD_x1 和 fFA_x1，这里赋值是为了避免重复计算以增加复杂度
[fMD_x1,fFA_x1] = f(P,x1); %重新计算新的x1对应的 fMD_x1 和 fFA_x1，以便下一轮循环使用
else %g(x1)>g(x2) 与上面情况相反
START_INT=x1;
x1 = x2;
x2=START_INT+tau*(END_INT-START_INT);
fMD_x1= fMD_x2;
fFA_x1 = fFA_x2;
[fMD_x2,fFA_x2] = f(P,x2);
end
k2=k2+1;
end
% 迭代结束后根据比较结果选择最终的解
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
strcmp函数：比较字符串（C语言中也有）
3.3 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.-->Ka随机且未知
%
% 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
% 输出是误检和虚警概率的误差底限 floor_MD 和 floor_FA。
% codebook size
M = 2^k;
%% 计算pbar
K_l = floor(E_Ka); K_u = ceil(E_Ka);
%E_Ka 的下取整和上取整，分别赋值给 K_l 和 K_u
while p_Ka(K_l-1) > 1e-12
K_l = K_l - 1;
end
K_l = max(K_l,0); %保证K_l大于等于0
while p_Ka(K_u+1) > 1e-12
K_u = K_u + 1;
end
%两个while循环的目的是寻找合适的区间[K_l,K_u]，使Ka落在该区间外的概率极小
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));
%% 初始化 floor_MD and floor_FA 为 pbar
floor_MD = pbar;
floor_FA = pbar;
%% Ka的求和
parfor Ka = K_l:K_u
%定义匿名函数来计算 Ka 和 KaEst(即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');%计算了从 KaEst_thres(Ka,K_l) 到 KaEst_thres(Ka,Ka-1) 的不完全伽马函数，'lower' 表示计算不完全伽马函数的下积分。
P_Ka_KaEst_list_b = gammainc(KaEst_thres(Ka,Ka+1:K_u),n,'upper');%计算了从 KaEst_thres(Ka,Ka+1) 到 KaEst_thres(Ka,K_u) 的不完全伽马函数，'upper' 表示计算不完全伽马函数的上积分。
P_Ka_Ka = 1 - max([P_Ka_KaEst_list_a P_Ka_KaEst_list_b]);
%计算了 P_Ka_Ka，它是上述两组概率的最大值与1的差。这个值在整个并行循环中用于计算误差底限
P_Ka_KaEst_list = [P_Ka_KaEst_list_a'; P_Ka_Ka; P_Ka_KaEst_list_b'];%将所有概率放入列表 P_Ka_KaEst_list 中
%加'是转置，所以列表矩阵矩阵 P_Ka_KaEst_list 的行数为 K_u - K_l + 1，列数为 3。其中，第一列是 P_Ka_KaEst_list_a 转置，第二列是 P_Ka_Ka，第三列是 P_Ka_KaEst_list_b 转置
%% 遍历Ka'(KaEst)
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;
% feval(p_Ka, Ka)：计算概率质量函数 p_Ka 在 Ka 处的值，表示活跃用户数为 Ka 的概率。
% max(Ka-KaEst_u,0)确保这个值不会为负
% P_Ka_KaEst 表示活跃用户数为 Ka 和 KaEst 之间的概率（从预先计算的概率矩阵 P_Ka_KaEst_list 中获取）
end
Mrx = (Ka + max(KaEst_l-Ka,0) - max(Ka-KaEst_u,0)); % 计算已解码的码字数量。
% max(KaEst_l-Ka,0) 这一项考虑了对于已经解码的活跃用户，可能会超过估计的范围 KaEst_l
% max(Ka-KaEst_u,0) 这一项考虑了对于未解码的活跃用户，可能会低于估计的范围 KaEst_u
if Mrx > 0
floor_FA = floor_FA + feval(p_Ka, Ka)*max(KaEst_l-Ka,0)*P_Ka_KaEst/Mrx;
% P_Ka_KaEst 表示在给定活跃用户数量 Ka 的情况下，估计用户数量在范围 [KaEst_l, KaEst_u] 内的概率。
end
end
end
floor_MD = min(floor_MD,1);
floor_FA = min(floor_FA,1);
%将误差底限限制在 [0, 1] 范围内
end
注: pbar 为论文[1]中的
；
公式为：
floor_MD 和 floor_FA 为论文[1]中的
；
feval函数：计算函数
3.4 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.
% Ka 服从泊松分布且未知，该函数用来计算低于给定的目标 misdetection（漏检）概率和 false alarm（误报）概率的最小所需 EbN0（每比特能量/单边噪声功率谱密度，以分贝为单位）。
% 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;
%% 调试模式
if DEBUG == 1
k = 128;
n = 19200;
epsilon_MD = .1; % Per-user error probability
epsilon_FA = .1; % Per-user error probability
E_Ka = 50;
rad_l = 0;
rad_u = 0;%采用 r=0 可获得最佳效果
end
%% 定义一个匿名函数来计算 Ka 的 PMF
% 允许用户使用不同的分布，只需替换这个函数定义
p_Ka = @(K) poisspdf(K,E_Ka);
%% 定义功率搜索范围
EbN0db_lowest = -2;%信噪比范围
EbN0db_highest = 15;
P_lowest = k*10^(EbN0db_lowest/10)/n;%对应的功率范围
P_highest = k*10^(EbN0db_highest/10)/n;
%% 寻找适当的解码半径（rad_l 和 rad_u）
%RCU_floor 函数用于计算误码率的底线，而目标是找到一个解码半径的范围，使得误码率底线低于给定的阈值 epsilon_MD 和 epsilon_FA
% 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
%这个额外的调整可能是为了确保系统在非常低的误码率条件下仍能正常工作，因为此时误码率目标相对于整个误码率范围来说可能较小
%% 定义匿名函数计算 RCU 界限
f_rcu = @(P,P1) RCU_KaRandomUnknown(P,P1,rad_l,rad_u,k,n,E_Ka,p_Ka);
%% 寻找最小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);
%% 储存结果
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
3.5 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.
%计算未知用户数量并遵循泊松分布系统的RCU
% 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;
%% 调试模式
if DEBUG == 1
k = 128;
n = 19200;
EbN0db = 4;
E_Ka_list = 50;
rad_l_list = [0];
rad_u_list = [0];
end
%% 码本
M = 2^k;
%% 符号功率预算
P_list = k.*10.^(EbN0db./10)./n;
%% 初始化
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));% 均为三维矩阵
%% 计算 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);%最内层循环计算不同条件下的 RCU
% 优化 P1，以最小化 max(p_MD,p_FA)
[p_MD_tmp,p_FA_tmp,P1_tmp] = golden_search_P1_MDFA(f_rcu,0,P,P/200,'max',[],[]);
% 最小化 p_MD and p_FA 的权重和, 用以下代码
% 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);
%% 储存结果
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
squeeze函数：删除长度为1的维度
3.6 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.
% 计算 Ka 随机且未知场景下的 PUPE 和 RCU
% 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
% 漏检和虚警的 RCU 边界
% codebook size
M = 2^k;
% number of samples for empirical evaluation of the CDF of It in qt
Niq = 1000;
%% 计算修正后的概率 ptilde
% 计算误差底限 [1, Th. 3]
[floor_MD,floor_FA] = RCU_floor_KaRandomUnknown(rad_l,rad_u,k,n,E_Ka,p_Ka);
%调用3.3的.m文件
% 寻找合适的 K_l 和 K_u ，使得 p_Ka 在区间 [K_l:K_u] 之外的概率很小。
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
%% 初始化 eps_MD and eps_FA 为 ptilde
eps_MD = ptilde;
eps_FA = ptilde;
%% 并行循环估计 Ka
parfor Ka = K_l:K_u
%% 计算 [1, Th. 2] 中的 xi_Ka_KaEst
% 计算 ML 估计器中 xi_Ka_KaEst 的参数
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); 见.m文件下
% 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_KaEst
xi_Ka_KaEst = P_Ka_KaEst_list(KaEst - K_l + 1);
% KaEst_l
KaEst_l = max(KaEst - rad_l,K_l);
% KaEst_u
KaEst_u = min(KaEst + rad_u,K_u);
% 创建两个集合 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);%t'下界
t1_u = max(KaEst_u-Ka,0) - max(KaEst_l-Ka,0);%t'上界
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
%% 计算 pt:
% 定义两个函数 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));
% 首先，考虑 t1 > 0 的情况
t1_vec2 = t1_vec;
t1_vec2(t1_vec2<0) = 0;
% 初始化 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]); %四维矩阵
% 预先计算一些量，后续使用
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;
% 求解三次方程来找到最大 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);% lambda 即为 λ
% 计算 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);
% 计算 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)]);
% 计算 ptt
ptt = min(exp(-n.*Ett_f),1);
% Now, 考虑 t1 = 0 的情况, 解一个二次方程找到最优的 λ
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);
% 结合 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
% 计算 pt ，即 ptt 数组每行的和
pt = min(sum(ptt,2),1);
%% 计算 t=1 时的 qt
qt = 1;
qtt = 1;
if t_vec(end) >= 1
t1_set = t1_vec2(2,:);
%由于计算复杂性问题，仅在 t = 1 且 Ka <= 150 的情况下计算 qt
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));% 生成大小为 (1, n) 的复高斯随机数 Zi
DefaultFA = sqrt(.5*max(KaEst_l-Ka,0)*P1)*(randn(1,n) + 1i*randn(1,n));%生成大小为 (1, n) 的默认 FA 部分的复高斯随机数 DefaultFA
codebook = sqrt(.5*P1)*(randn(min(Ka,KaEst_u),n) + 1i*randn(min(Ka,KaEst_u),n));%生成大小为 (min(Ka, KaEst_u), n) 的码本的复高斯随机数 codebook
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);%取 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
% 取 pt and qt 的最小值
pt(2) = min(pt(2),qt);
ptt(2,:) = min(ptt(2,:),qtt);
end
% 分别取 pt, qt 与 xi(Ka,Ka') 的最小值
pt = min(pt,xi_Ka_KaEst);
ptt = min(ptt,xi_Ka_KaEst);
%% 为给定的 P' < P 计算 epsilon_MD 和 epsilon_FA
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)); % 已解码码字的数量
Mrx(Mrx == 0) = inf; %将 Mrx 中的零元素替换为无穷大
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);%得到epsilon_MD 和 epsilon_FA
end
注：修正后的概率 ptilde即为论文[1]中的
；
公式为：
注： xi_Ka_KaEst 即为论文[1]中的
;
KaEst_thres即为论文[1]中的
；
公式为：
对于基于能量的估计器(即文章中的 For energy-based estimator of Ka )
注： KaEst_l 和 KaEst_u 分别是论文[1]中的
；
预计算的量：
注： ptt 表示论文[1]中的
；
公式为：
reshape函数：重构数组（与repmat函数不同）
注： qt 和 qtt 即为论文[1]中的
；
；
公式为：
注： .’的含义是非共轭转置
四、SA-MPR
4.1 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.
%计算在 Ka 服从泊松分布且已知的条件下，通过使用具有多包接收（MPR）的时隙ALOHA（SA）实现的PUPE（用户包错误概率），并找到使其低于给定阈值epsilon所需的最小EbN0
% 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 指示是否使用 slot-index coding 的标志，为1表示使用，为0表示不使用
%
% OUTPUT
% data : store the system parameters and the minimal required EbN0
addpath ../RCU_KaKnown %将路径 "../RCU_KaKnown" 添加到 MATLAB 的搜索路径中。这是为了确保 MATLAB 能够找到位于这个路径下的函数和文件。
tic
DEBUG = 0;
if DEBUG == 1
k = 128;
n = 19200;
epsilon = .1;
E_Ka = 50;
end
p_Ka = @(K) poisspdf(K,E_Ka);
%% 搜索的功率范围
EbN0db_lowest = -1;
EbN0db_highest = 15;
P_lowest = k*10^(EbN0db_lowest/10)/n;
P_highest = k*10^(EbN0db_highest/10)/n;
%% 计算 RCU （与之前不同，考虑是否使用SlotIdxCoding）
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 %区别在于k是否要减去floor(log2(L))，其中 PMF_Ksa函数 定义在下面
%% 寻找所需最小 EbN0
% 通过 P1 和 L 对 RCU 进行优化
[eps_RCU,P_RCU,P1,L] = binary_search_P_SAMPR_KaKnown(f_rcu, P_lowest, P_highest,...
epsilon,epsilon/100,E_Ka);%其中 binary_search_P_SAMPR_KaKnown函数定义在下面
EbN0db_RCU = 10*log10(n*P_RCU/k);
%% 储存结果
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)
% 计算每个时隙中 Ka 的 PMF
% Ksa 是待计算PMF的活跃用户数量的向量
PMF_Ksa = zeros(size(Ksa));
K_u = E_Ka;
while p_Ka(K_u+1) > eps
K_u = K_u + 1;
end%设置一个上限 K_u，使得概率质量函数 p_Ka 在 K_u+1 处的值小于阈值 eps ，eps 表示接近于零的一个小数。这个循环用于找到 p_Ka 非零的最大 K_u
K_u = max(K_u,10000);%最小也要10000
PMF_Ksa_scalar = @(T) sum(p_Ka(0:K_u).*binopdf(T,0:K_u,1/L));
%计算对应于 T 的 binopdf 和 p_Ka 的乘积之和。这表示计算了在每个时隙中有 T 个活跃用户的概率。
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;
k1=0;
[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);
% golden_search_P1_SAMPR_KaKnown函数定义在下面，在给定范围内搜索参数 P1 和 nSlot，使得目标函数 f 在 P1 和 nSlot 取值处的结果最小
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)
% 在给定 P 的情况下，最小化 P1 上的 min_L f(P,P1,L)
P = END_INT;
iter= 20;
tau=double((sqrt(5)-1)/2); % 0.618
k2=0;
x1=START_INT+(1-tau)*(END_INT-START_INT);
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);
% golden_search_nSlot_SAMPR_KaKnown函数定义在下面
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)
% 在给定 P 和 P1 的情况下，找到最小 f(P,P1,L) 和 L
iter= 20;
tau=double((sqrt(5)-1)/2); % 0.618
k2=0;
x1=max(round(START_INT+(1-tau)*(END_INT-START_INT)),0);
x2=round(START_INT+tau*(END_INT-START_INT));
f_x1 = f(P,P1,x1);
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
注： Slot-Index Coding（时隙索引编码）是一种多用户通信中的编码技术，特别用于解决随机接入通信系统中的干扰问题。在随机接入通信中，多个用户在同一时隙内传输信息，可能导致碰撞和干扰。时隙索引编码通过在传输时指定时隙索引，以降低干扰水平。
在时隙索引编码中，每个用户都被分配一个唯一的时隙索引，并且通过选择适当的编码方式，用户可以在同一时隙内传输信息而不会相互干扰。这有助于提高通信系统的性能，特别是在高负载和高干扰的环境中。
4.2 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
%对于 Ka 服从泊松分布且未知的情况
% 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 %与泊松已知的情况区别在 DecRad 上，DecRad 是解码半径
%% 寻找最小所需的 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); % binary_search_P_SAMPR 函数定义见下
EbN0db_RCU = 10*log10(n*P_RCU/k);
%% 储存结果
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)
% 计算每个时隙内 Ka 的 PMF
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;
k1=0;
[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);% golden_search_P1_SAMPR 函数定义见下
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)
% 优化时隙数
iter= 20;
tau=double((sqrt(5)-1)/2); %0.618
k2=0;
x1=max(round(START_INT+(1-tau)*(END_INT-START_INT)),0);
x2=round(START_INT+tau*(END_INT-START_INT));
%round函数见下
[fMD_x1,fFA_x1,DecRad_x1] = linear_search_DecRad_SAMPR(f,P,P1,x1,weight_MD,weight_FA); % computing values in x points linear_search_DecRad_SAMPR 函数定义见下
[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 %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)
% 优化解码半径
maxrad = 5; % 从0-4中选择解码半径
fMD = zeros(maxrad,1);
fFA = zeros(maxrad,1);
for ii = 1:maxrad
[fMD_tmp, fFA_tmp] = f(P,P1,x,ii-1);%存储不同解码半径下的 fMD 和 fFA
fMD(ii) = fMD_tmp;
fFA(ii) = fFA_tmp;
end
[~,idxMin] = min(fMD*weight_MD + fFA*weight_FA);
% ~ 用于忽略 min 函数的输出，只关心最小值的下标 idxMin ，而不需要获取最小值本身
fMD = fMD(idxMin);
fFA = fFA(idxMin);
iiset = 0:maxrad-1;
DecRad = iiset(idxMin); % 解码半径为0到4
% DecRad = 0;
% [fMD, fFA] = f(P,P1,x,DecRad);
% 意思是 也可以直接将 DecRad 设为 0，然后调用 f 函数获取对应的 fMD 和 fFA
end
round函数：舍入到最近的小数或者整数
五、TIN
5.1 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
% 在把 干扰是为噪声(TIN) ，且Ka服从泊松，未知的情况下，计算达到 epsilon_MD, epsilon_FA 所需的最小EbN0
% 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;
%% 定义匿名函数来计算 RCUs-TIN 界限
f_TIN = @(P,P1,s) RCUs_TIN(P,P1,s,k,n,E_Ka,p_Ka,normalApprox);% RCUs_TIN 函数定义见下
%% 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);% binary_search_TIN 函数定义见下
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 这里的n和上面不一样，这里是块长度
% P = transmit power (no dB!) 传输功率，单位不是dB！
% R = Values of Rate ln(M)/n at which the bound will be computed 以 ln(M)/n 的比值来计算界限
% Ka = Number of active users
% s = parameter to optimize in the RCUs (s>0) 在 RCUs 中进行优化的参数，s=1时使用正态近似
% 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)); % nan函数 见.m代码下
parfor ii = 1:length(Ka_vec) %并行循环
Ka = Ka_vec(ii);%获取当前循环迭代中的活跃用户数（Ka）
snr = P/(1+(Ka-1)*P);%计算信噪比
x = sqrt(0.5*snr)*(randn(n,N) + 1i*randn(n,N)); % 生成一个信号由其中一个用户发送。
z = sqrt(0.5)*(randn(n,N) + 1i*randn(n,N)); % 生成噪声和干扰
y = x + z; % 接收信号
i_s = -s*vecnorm(y-x).^2 + s*vecnorm(y).^2/(1+s*snr) + n*log(1+s*snr); % generalized information density
% vecnorm函数 定义见.m代码下
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 ptilde
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)); %论文4页公式(8)，见.m代码下
%% Initialize eps_MD and eps_FA to be ptilde
eps_MD = ptilde;
eps_FA = ptilde;
%% Compute the xi_Ka_KaEst 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 %s等于1时，使用正态近似
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)));% qfunc函数和eta_s_vec公式 见.m代码下
else
eta_s_vec = eta_s(n,P1,R,(K_l:K_u)',s,N);%s不等于1时
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;
k1=0;
[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);% golden_search_P1_TIN函数 定义见下
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;
tau=double((sqrt(5)-1)/2); % g0.618
k2=0;
x1=START_INT+(1-tau)*(END_INT-START_INT); % computing x values
x2=START_INT+tau*(END_INT-START_INT);
s_start = 0;% s 的范围
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);
% golden_search_s_TIN函数 定义见下，用来搜索最优的参数 s_opt1 和 s_opt2
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;
tau=double((sqrt(5)-1)/2); % 0.618
k2=0;
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
generalized information density：广义信息密度，计算公式：
信噪比（SNR）：计算公式如下
NaN函数：创建所有值均为 NaN 的数组
vecnorm函数：向量范数（默认p=2）
注: ptilde公式为
注： \eta_s即为论文[1]中的
；
公式为：
s等于1时，
s不等于1时，
注： C_TIN 和 V_TIN 的计算公式为
qfunc函数： Q函数
注： eps_MD公式为
eps_FA公式为
六、论文中仿真图
6.1 Eb/N0—E[Ka]
为了达到 max{PMD, PFA} ≤ 10−1，所需的 Eb/N0 作为 E[Ka] 的函数，其中 k = 128 bits，n = 19200 信道，Ka ∼ Pois(E[Ka])。实线代表 Ka 未知的方案/界限；虚线代表 Ka 已知的方案/界限
6.2 MD与FA概率—Eb/N0
MD和FA概率的界限作为Eb/N0的函数，其中k = 128位，n = 19200 信道
6.3 Eb/N0—E[Ka]
为了达到 {PMD ≤ 10−1 , PFA ≤ 10−3 }，{PMD ≤ 10−3 , PFA≤ 10−1}，或者 max{PMD, PFA} ≤ 10−3 作为E[Ka]的函数，需要的Eb/N0，其中k = 128比特，n = 19200 信道，且Ka ∼ Pois(E[Ka])。
6.4 MD与FA概率—Eb/N0
MD和FA概率的界限作为Eb/N0的函数，对于k = 128位和n = 19200信道，选择与6.3相对应的rℓ和ru。实线代表εMD；虚线代表εFA。
其余图暂不展示
七、论文复现图片
7.1 EbN0_KaPoissonKnown(100,15000,0.1,2)
程序:
%E7%A8%8B%E5%BA%8F%E5%9B%BE.png)
过程:
%E8%BF%87%E7%A8%8B%E5%8F%98%E9%87%8F.png)
结果:
%E7%BB%93%E6%9E%9C%E5%9B%BE.png)
7.2 RCU_KaPoissonKnown(128,19200, 50, 0)
过程:
%E8%BF%87%E7%A8%8B%E5%8F%98%E9%87%8F.png)
程序和结果:
%E7%BB%93%E6%9E%9C%E5%9B%BE.png)
7.3 RCU_KaFixedKnown(100,15000,100,2)
程序和结果：
.png)
7.4 EbN0_KaPoissonUnknown
跑不出来
7.5 RCU_KaPoissonUnknown
跑不出来
7.6 EbN0_SAMPR_KaPoissonKnown( 128, 19200, 0.1, 50, 1)
过程：
%E8%BF%87%E7%A8%8B%E5%8F%98%E9%87%8F.png)
程序和结果：
%E7%BB%93%E6%9E%9C%E5%9B%BE.png)
7.7 EbN0_SAMPR_KaPoissonUnknown( 128, 19200, 0.1, 0.1, 50, 0)
过程：
%E8%BF%87%E7%A8%8B%E5%8F%98%E9%87%8F.png)
结果：
%E7%BB%93%E6%9E%9C%E5%9B%BE.png)
7.8 EbN0_TIN_KaPoissonUnknown(128, 19200,0.1,0.1, 100,1)
过程：
%E8%BF%87%E7%A8%8B%E5%8F%98%E9%87%8F.png)
程序和结果：
%E7%BB%93%E6%9E%9C%E5%9B%BE.png)
博客来自：祝小火💗，后续会持续改动~~
【推荐】国内首个AI IDE，深度理解中文开发场景，立即下载体验Trae
【推荐】编程新体验，更懂你的AI，立即体验豆包MarsCode编程助手
【推荐】抖音旗下AI助手豆包，你的智能百科全书，全免费不限次数
【推荐】轻量又高性能的 SSH 工具 IShell：AI 加持，快人一步
· 全程不用写代码，我用AI程序员写了一个飞机大战
· DeepSeek 开源周回顾「GitHub 热点速览」
· MongoDB 8.0这个新功能碉堡了，比商业数据库还牛
· 记一次.NET内存居高不下排查解决与启示
· 白话解读 Dapr 1.15：你的「微服务管家」又秀新绝活了