标签：Aided decoded Codes end Unsourced current num round roots

一、摘要

随记仅用于个人对论文的分析、初步复现。

1.1 文件夹介绍

随机包含了一篇论文的仿真结果的源代码，该论文的标题是 "CRC-aided Spare Regression Codes for Unsourced Random Access"。

源代码CRC-aided_SPARCs_for_URA-main，一共包括三个文件夹：

"CRC-BMST codes for stitching"，

**"CRC-aided SPARCs for URA" **，

"parameter configurations"。

关于 "CRC-BMST codes for stitching" 文件夹，它包括：

1.CRC-BMST codes文件夹，其中包含了我们提出的 CRC-BMST 码的 MATLAB 源代码的实现；

2.Alex_Fengler_tree_codes文件夹，其中包含了实现原始树码的 MATLAB 源代码（由 Alex Fengler 博士撰写）。

关于 "CRC-aided SPARCs for URA" 文件夹，它包括用于通过我们提出的级联编码方案模拟整个无源随机访问信道模型的 MATLAB 源代码。

关于 "parameter configurations" 文件夹，它包含了三个存储论文仿真结果中使用的具体参数的 '.mat' 文件。

1.2 m文件介绍

CRC-BMST codes for stitching文件夹中的m文件：

CRC-BMST codes子文件夹：

1.CRC_check.m

2.CRC_encoding.m

3.outer_CRC_checks_BMST_efficient.m

4.outer_encoding_BMST.m

5.UMAC_lossless_stitching_CRC_BMST.m

Alex_Fengler_tree_codes子文件夹：

1.outerDecoder.m

2.outerEncoder.m

3.UMAC_lossless_stitching_tree_codes.m

CRC-aided SPARCs for URA文件夹中的m文件：

1.CRC_check.m

2.CRC_encoding.m

3.outer_CRC_checks_BMST_efficient.m

4.outer_encoding_BMST.m

5.Outer_UMAC_stitching_CRC_BMST.m

6.AMP_SPARCs_UMAC.m

7.MAP_AMP_Hybrid_SPARCs_UMAC.m

8.MAP_AMP_hybrid_UMAC_CRC_BMST.m

9.MAP_AMP_hybrid_UMAC_CRC_BMST_extra_iteration.m

10.FWHT_user_subsampling.m

另外，该文件夹还包含一个.c 文件和一个mexw64文件

其中.mexw64 文件是 MATLAB 编译的二进制 MEX 文件的一种类型。MEX 文件是用 C、C++ 或 Fortran 编写的 MATLAB 可执行函数，通过编译生成的二进制文件。.mexw64 文件是在 Windows 平台上生成的 MEX 文件，并且是 64 位 Windows 系统的版本。

11.fwht_user.c

12.fwht_user.mexw64

1.3 mat文件介绍

parameter configurations文件夹中的mat文件：

1.configurations_UMAC_J_15_L_11_B_75 (Fig_3).mat

2.configurations_UMAC_J_14_L_7_B_50 (Fig_4).mat

3.configurations_UMAC_J_15_L_11_B_75 (Fig_5).mat

二、CRC-BMST codes for stitching

2.1 CRC-BMST codes

实现框图：

2.1.1 CRC_check.m

function check = CRC_check(r_c, g)% 该函数的目的是通过 CRC 校验来确定接收到的多项式是否包含错误
check = 0;
[~, r] = gfdeconv(fliplr(r_c), fliplr(g));% r_c 是接收到的多项式，g 是生成多项式。gfdeconv 函数返回商和余数，但在此处只取了余数，即 r
if (sum(r)==0)
    check = 1;
end % 检查余数 r 的各项之和是否为零。如果为零，说明整个除法过程没有余数，即接收到的多项式是没有错误的。此时，将 check 设置为 1，表示通过了 CRC 校验

end

注：1、gfdeconv函数：用于对两个有限域（GF）多项式进行除法运算

[q, r] = gfdeconv(a, b, m)

其中，a和b是GF多项式的系数向量，m是GF的阶数（GF(2^m)）。

函数返回商多项式q和余数多项式r，使得a = conv(b, q) + r，其中conv表示多项式的卷积运算。m默认为1，GF(2).

2、fliplr 函数:将数组从左向右翻转，因为 gfdeconv 函数要求输入的多项式系数是从高次到低次排列的

3、关于CRC码：是采用二进制除法(没有进位，用XOR来代替减法)相除得到的余数，在做除法之前，要在信息数据之后先加上r个0.

2.1.2 outer_CRC_checks_BMST_efficient.m

function [current_blocks, check] = outer_CRC_checks_BMST_efficient(previous_blocks, previous_superimposed_sum, current_path, current_block, num_protect_section, m, memory, r, CRC_poly)
% 用于进行 CRC（Cyclic Redundancy Check）检查的程序，其中使用了 BMST（Block Markov Superposition Transmission）编码的结果
% store previous original_block_bits so that the same info doesn't need to compute twice.
% start from the second layer.

% current_blocks：存储当前阶段（或轮次）的恢复信息比特。
% check：指示当前比特是否满足 CRC 条件的指示器。

% previous_blocks：所有先前恢复的信息块比特。
% previous_superimposed_sum：先前计算的比特级别的累积和,与current_path结合可获取当前阶段的累积和。
% current_path：包含所有节点的当前路径。
% current_block：当前阶段或轮次。
% num_protect_section：受保护部分的数量。
% m：每个块的比特数，即 log2(M)。
% memory：在 BMST 编码中使用的编码记忆参数。
% r：指示每轮的冗余比特数量的向量。
% CRC_poly：选择的 CRC 生成多项式。
% 这些参数主要用于BMST，可不详细研究

    check = 0;
    original_block_bits = zeros(current_block, m);% 初始化矩阵，其中包含了先前阶段所有块的信息比特
    original_block_bits(1:current_block-1, :) = previous_blocks;% 将先前已经恢复的块的信息比特复制到 original_block_bits 中。    
    current_block_bits = de2bi(current_path(current_block)-1, m);%将当前阶段的信息块下标转换为二进制，并存储在 current_block_bits 中。
    original_block_bits(current_block, :) = mod(current_block_bits-previous_superimposed_sum, 2); % 计算并存储当前阶段的信息比特。这里，对 current_block_bits 减去 previous_superimposed_sum，然后取模 2 (实际上是异或运算（XOR）)

    current_path_info_bits = [];% 用于存储当前路径的信息比特
    start_cc = 0;% 用于追踪当前路径信息比特的起始位置

    for cc = current_block - num_protect_section+1 : current_block-1  % 循环遍历当前路径中受保护的部分
      end_cc = start_cc + m - r(cc);% 计算当前路径信息比特的结束位置,r(cc)表示本轮的冗余。
      current_path_info_bits(start_cc+1: end_cc) = original_block_bits(cc, 1:m-r(cc));% 将先前阶段中不包含冗余比特的信息比特复制到 current_path_info_bits 中，形成当前路径的信息比特。
      start_cc = end_cc;% 更新 start_cc 为当前路径信息比特的结束位置，以便下一次循环计算。
    end

    current_blocks = original_block_bits;

    check_bits = [current_path_info_bits, original_block_bits(current_block, :)];% 将当前路径中的信息位current_path_info_bits和当前块的信息位original_block_bits(current_block, :)合并成一个检查位序列check_bits
    check = CRC_check(check_bits, CRC_poly);
    % 调用CRC_check函数进行CRC检查，将check_bits和CRC生成多项式CRC_poly传递给该函数，检查是否通过CRC验证。如果通过，check被设置为1，表示检查通过；否则，check为0，表示检查未通过。


end

%思路就是先对接收块实现BMST，最后用CRC来检查是否出错

注：1、概念

BMST 编码（Block Markov Superposition Transmission）： BMST 是一种用于通信系统的编码技术。在 BMST 编码中，信息比特被分块处理，每个块的编码依赖于前一块的编码结果，这种依赖关系是通过马尔可夫链实现的。BMST 编码可以用于提高通信系统的性能和可靠性。
CRC 检查（Cyclic Redundancy Check）： CRC 是一种纠错校验方法，通常用于检测或纠正数据传输中的错误。CRC 使用多项式除法来生成冗余校验码，然后将这些校验码附加到原始数据中。在接收端，通过再次进行多项式除法，可以检测或纠正数据中的错误。

2、de2bi函数：de2bi 函数用于将这个整数转换为二进制表示，并 m 表示每个二进制表示的比特数。

例如，如果 m 是 4，那么 de2bi(current_path(current_block)-1, m) 将把 current_path(current_block)-1 转换为一个包含 4 位的二进制向量。这个二进制向量表示了信息块的信息比特。

2.1.3 CRC_encoding.m

function c = CRC_encoding(m, g)
% m ：信息向量
% g ：生成多项式
% c : CRC编码的结果
% 注意：由于gfdeconv中的多项式次数是按升序排列的，因此需要执行fliplr操作。

c = zeros(1, length(g)-1);
m_r = [m zeros(1, length(g)-1)];% 在信息向量 m 的末尾添加零(后面留的是冗余位)
[~, r] = gfdeconv(fliplr(m_r), fliplr(g));% m_r 除以 g 得到的CRC码的长度会比 g 的长度大，因为需要加冗余位，这里只保留余数 r (r也是向量)
c(1:length(r)) = r;% 将余数 r 的内容复制到CRC编码结果 c 中
c = fliplr(c);% 由于之前进行了fliplr操作，这里再次使用fliplr将结果翻转，以得到正确的CRC编码结果。

end

2.1.4 outer_encoding_BMST.m

function outer_encoded = outer_encoding_BMST(message, memory, m, num_round, r, protect_sections, num_info_bits_round, orderings)
  %  CRC-BMST 编码过程
  %  只有信息比特（message bits）受到CRC码的保护，而不是全部比特。

  % message: 要进行编码的消息比特;
  % memory: BNST中的编码内存;
  % m: 每块的比特数;
  % num_round: 每个用户的原始消息被分割的次数;
  % r: 一个向量，指示每一轮的冗余比特数;
  % protect_sections: 前面受保护的块的数量;
  % num_info_bits_round: 一个向量，指示每一轮的信息比特数;
  % orderings: CRC-BMST编码过程中使用的交织器;


% 不同CRC位下的CRC生成多项式(r位CRC码最多可保护2^r-1-r位信息位)
poly{4} = [1 0 0 1 1]; % 信息位为2^4-1-4=11

%poly{5} = [1 0 0 1 0 1]; % 信息位为 26
poly{5} = [1 0 1 0 1 1]; % 信息位为 10

%poly{6} = [1 0 0 0 0 1 1]; % 信息位为 57
poly{6} = [1 0 0 0 1 1 1]; % 信息位为 25

poly{7} = [1 0 1 1 0 1 1 1]; % 信息位为 56

%poly{8} = [1 0 0 1 0 1 1 1 1]; % 信息位为 119
%poly{8} = [1 0 0 0 0 0 1 1 1]; % 信息位为 119
poly{8} = [1 1 1 0 1 0 1 1 1]; % 信息位为 9

%poly{9} = [1 0 1 0 0 1 0 1 1 1]; % 信息位为 246
%poly{9} = [1 0 1 1 1 1 1 0 1 1]; % 信息位为 246
%poly{9} = [1 1 0 0 0 0 1 0 1 1]; % 信息位为 13
poly{9} = [1 0 0 1 1 1 1 0 0 1]; % 信息位为 8

%poly{10} = [1 1 0 0 0 1 1 0 0 1 1]; % 信息位为 501
%poly{10} = [1 0 1 0 1 1 1 0 0 1 1]; % 信息位为 21
%poly{10} = [1 0 1 0 0 0 1 1 1 0 1]; % 信息位为 12
poly{10} = [1 0 1 0 0 1 1 0 1 1 1]; % 信息位为 5

% poly{11} = [1 0 0 1 1 1 1 0 1 0 1 1]; % 信息位为 4
poly{11} = [1 0 1 0 1 1 1 0 0 0 1 1]; % 信息位为 12

poly{12} = [1 0 1 0 0 1 0 0 1 1 1 1 1]; % 信息位为 11

poly{13} = [1 0 0 0 0 1 0 1 1 0 1 1 1 1]; % 信息位为 11

poly{14} = [1 0 0 0 1 1 0 1 1 1 0 0 0 1 1]; % 信息位为 11


  outer_encoded(1:m) = message(1:m);%  message 的前 m 位直接赋值给 outer_encoded，因为信息位不需要进行 CRC-BMST 编码。
  message_last_end = m;% 表示上一轮编码结束的位置，初始值为 m
  original_encoded_block = zeros(num_round, m);% 用于存储每一轮的编码块，初始值全部为零
  original_encoded_block(1, :) = message(1:m);% message 的前 m 位设为第一轮的原始信息块 

  for current_round = 2 : num_round% 从第二轮开始进行循环，生成每一轮的编码块
    protected_info_bits = [];% 用于存储受 CRC 保护的信息位
    num_protect_section = protect_sections(current_round);% 获取当前轮的保护区段数
    encoded_last_end = (current_round-1)*m;% 计算上一轮编码结束的位置(一轮m位)
    info_current = message(message_last_end+1 : message_last_end + m -r(current_round));% 从 message 中提取当前轮需要编码的信息位
    start_ = 0;
    if (current_round > num_protect_section)% 如果当前轮大于保护区段数
      start_ = num_info_bits_round(current_round-num_protect_section);% 更新 start_ 为前面所有轮的信息位总数
    end
    protected_info_bits =message(start_+1 : num_info_bits_round(current_round));% message中提取出当前轮需要保护的信息位。
    outer_encoded_current = [protected_info_bits, CRC_encoding(protected_info_bits, poly{r(current_round)})];% 保护的信息位和CRC编码后的冗余位拼接在一起，得到当前轮的编码块

    len_outer_current = length(outer_encoded_current);% 计算当前轮编码后的长度
    start_block_BMST = max(current_round-memory, 1);% 计算BMST编码的起始位置，确保不小于1
    outer_encoded_block = outer_encoded_current(len_outer_current-m+1:len_outer_current);% 从当前编码后的序列中提取出当前轮的编码块
    original_encoded_block(current_round, :) = outer_encoded_block;% 当前轮的编码块存储到 original_encoded_block 中
    num_m = 0;% 用于记录BMST编码过程中的迭代次数
    for sb = current_round-1 : -1 : start_block_BMST% 从当前轮的前一轮开始，逐轮进行BMST编码，直到start_block_BMST轮
      num_m = num_m + 1;
      permuting_block = original_encoded_block(sb, :);% 取出原始编码块，准备用于BMST编码
      outer_encoded_block = mod(outer_encoded_block + permuting_block(orderings(num_m, :)),2);% 将当前BMST迭代的结果加到已编码块上
    end
    outer_encoded(encoded_last_end+1:encoded_last_end+m) = outer_encoded_block;% 将当前轮的已编码块存储到总的已编码序列中
    message_last_end = message_last_end + m - r(current_round);% 更新当前消息的结束位置
  end

end

注：

2、poly函数：

2.1.5 UMAC_lossless_stitching_CRC_BMST.m

function [PUPE, average_num_checks] = UMAC_lossless_stitching_CRC_BMST(num_active, num_round, num_info_bits, M, r, protect_sections, num_iterations, memory)
  % use block Markov superposition transmission to (lossless) stitch chunks

  % PUPE:  per-user probablity of error;
  % average_num_checks: the average number of checks computated during stitch;

  % num_active: the number of active users;
  % num_round: the original messages transmitted by each user are splitted into such rounds;
  % M : the section size, i.e., M = 2^m where m is the number of bits each round;
  % r: a vector indicating the number of redundant bits at each round;
  % protect_sections: the number of preceding sections protected at current round;
  % num_iterations: the number of simulation runs;
  % memory: the encoding memory of BMST.

  
  % the CRC generator polynomials with different number of CRC bits
poly{4} = [1 0 0 1 1]; % the number of info bits up to length 11

% poly{5} = [1 0 0 1 0 1]; % up to length 26
poly{5} = [1 0 1 0 1 1]; % up to length 10

% poly{6} = [1 0 0 0 0 1 1]; % up to length 57
poly{6} = [1 0 0 0 1 1 1]; % up to length 25

poly{7} = [1 0 1 1 0 1 1 1]; % up to length 56

% poly{8} = [1 0 0 1 0 1 1 1 1]; % up to length 119
% poly{8} = [1 0 0 0 0 0 1 1 1]; % up to length 119
poly{8} = [1 1 1 0 1 0 1 1 1]; % up to length 9

% poly{9} = [1 0 1 0 0 1 0 1 1 1]; % up to length 246
% poly{9} = [1 0 1 1 1 1 1 0 1 1]; % up to length 246
% poly{9} = [1 1 0 0 0 0 1 0 1 1]; % up to length 13
poly{9} = [1 0 0 1 1 1 1 0 0 1]; % up to length 8

% poly{10} = [1 1 0 0 0 1 1 0 0 1 1]; % up to length 501
% poly{10} = [1 0 1 0 1 1 1 0 0 1 1]; % up to length 21
% poly{10} = [1 0 1 0 0 0 1 1 1 0 1]; % up to length 12
% poly{10} = [1 0 1 0 0 1 1 0 1 1 1]; % up to length 5

% poly{11} = [1 0 0 1 1 1 1 0 1 0 1 1]; % up to length 4
poly{11} = [1 0 1 0 1 1 1 0 0 0 1 1]; % up to length 12

poly{12} = [1 0 1 0 0 1 0 0 1 1 1 1 1]; % up to length 11

poly{13} = [1 0 0 0 0 1 0 1 1 0 1 1 1 1]; % up to length 11

poly{14} = [1 0 0 0 1 1 0 1 1 1 0 0 0 1 1]; % up to length 11

L = num_active;
m = log2(M);


% randomly generate interleavers used for BMST
orderings = zeros(memory, m);
inverse_orderings = zeros(memory, m);
for mmr = 1 : memory
  perm_current = randperm(m);
  inverse_perm_current(perm_current) = 1 : m;
  orderings(mmr, :) = perm_current;
  inverse_orderings(mmr, :) = inverse_perm_current;
end



num_info_bits_round(1) = m;

for current_round = 2 : num_round
  num_info_bits_round(current_round) = num_info_bits_round(current_round-1) + m - r(current_round);
end

num_remains = zeros(1, num_iterations);
num_candidates = num_active;
average_checks = zeros(1, num_iterations);
for iteration = 1 : num_iterations
   rng('shuffle');
   num_errors_current = 0;
   outer_encoded = zeros(L, m*num_round);
   outer_encoded_pos = zeros(num_round, L);
   message = zeros(L, num_info_bits);

   y = cell(1, num_round);

   % outer encoding via CRC codes at each round
   for l = 1 : L
     message(l, :) = randi(2, 1, num_info_bits)-1;
     if (sum(message(l, 1:m))==0)
       message(l, m) = 1;
     end
     outer_encoded(l, :) = outer_encoding_BMST(message(l, :), memory, m, num_round, r, protect_sections, num_info_bits_round, orderings);
     outer_encoded_pos(:, l) = bi2de(reshape(outer_encoded(l, :), [m, num_round])')+1;
   end

   % (lossless) stitching procedure
   decoded_CRC_encoded_pos = outer_encoded_pos;
   
   % count duplicate roots
   decoded_roots_msg= decoded_CRC_encoded_pos(1, :);

   num_decoded_duplicate_roots = 0;
   duplicate_decoded_roots = {};
   [~, unique_roots_idx] = unique(decoded_roots_msg, 'stable');
   roots_duplicate_idx = setdiff(1:num_candidates, unique_roots_idx);
   
   while (~isempty(roots_duplicate_idx))
       num_decoded_duplicate_roots = num_decoded_duplicate_roots + 1;
       current_root_idx = roots_duplicate_idx(1);
       current_root_msg = decoded_roots_msg(current_root_idx);
       duplicate_decoded_roots_current = find(decoded_roots_msg==current_root_msg);
       duplicate_decoded_roots{num_decoded_duplicate_roots} = duplicate_decoded_roots_current;
       roots_duplicate_idx = setdiff(roots_duplicate_idx, duplicate_decoded_roots_current);
   end
   
   
   
   num_survival_paths = zeros(num_candidates, num_round);
   survival_current_blocks = cell(num_candidates, num_candidates);

   survival_paths = cell(num_candidates, num_round);  % each cell represents a matrix containing all survival paths.
   
   for current_user = 1 : num_candidates
     survival_paths{current_user, 1} = decoded_CRC_encoded_pos(1, current_user);
     num_survival_paths(current_user, 1) = 1;
     survival_current_blocks{current_user, 1} = de2bi(decoded_CRC_encoded_pos(1, current_user)-1, m);
   end


   for current_round = 2 : num_round
       survival_blocks_prev = cell(num_candidates, num_candidates);
       survival_blocks_prev = survival_current_blocks;
       survival_current_blocks = cell(num_candidates, num_candidates);
       for current_root = 1 : num_candidates
           num_paths_current = 0;
           num_paths_prev = num_survival_paths(current_root, current_round-1);
           survival_paths_prev = survival_paths{current_root, current_round-1};
           survival_path_current = [];
           unique_msg_current = unique(decoded_CRC_encoded_pos(current_round, :), "stable");
           num_protect_section = protect_sections(current_round);
           for current_path_num = 1 : num_paths_prev
               survival_path_prev = survival_paths_prev(current_path_num, :);
               previous_blocks = survival_blocks_prev{current_root, current_path_num};
               superimposed_sum = zeros(1, m);
               num_m = 0;
               start_block_BMST_ = max(current_round-memory, 1);
               for prev_batch = current_round-1 : -1 : start_block_BMST_
                   num_m = num_m + 1;
                   superimposed_sum = superimposed_sum + previous_blocks(prev_batch, orderings(num_m,:));
               end
               
               previous_superimposed_sum = mod(superimposed_sum, 2);
               for current_pos_considered = unique_msg_current
                   CRC_poly = poly{r(current_round)};
                   current_path = [survival_path_prev, current_pos_considered];
                   current_blocks = [];
                   [current_blocks,check] = outer_CRC_checks_BMST_efficient(previous_blocks, previous_superimposed_sum, current_path, current_round, num_protect_section, m, memory, r, CRC_poly);
                   if (check)
                       num_paths_current = num_paths_current + 1;
                       survival_path_current(num_paths_current, 1:current_round) = [survival_path_prev, current_pos_considered];
                       survival_current_blocks{current_root, num_paths_current} = current_blocks;
                   end
               end
           end
           num_survival_paths(current_root, current_round) = num_paths_current;
           survival_paths{current_root, current_round} = survival_path_current;
       end
   end
   
   
   % peeling decoding procedure in the sense that peel off nodes included in the unique path
   remaining_message = cell(1, num_round);
   for current_round = 1 : num_round
     remaining_message{current_round} = decoded_CRC_encoded_pos(current_round, :);
   end

   num_remainings_last = num_candidates;
   remaining_roots_elimination = 1 : num_candidates;
   average_checks(iteration) = sum(sum(num_survival_paths(:,1:num_round-1), 2));
   final_num_survival_paths = num_survival_paths(:, num_round);
   num_detected_error_trees = sum(final_num_survival_paths==0);
   remaining_roots_elimination(final_num_survival_paths==0) = []; % excluding the detected failure trees.
   num_decoded_users = 0;

while (1)     
        
   while (1)
        num_remainings = 0;
        remaining_roots = [];
        current_remaining_roots_elimination = remaining_roots_elimination;

        % eliminating the corrected messages
        for current_root  = remaining_roots_elimination
          if (num_survival_paths(current_root, num_round)==1)
            flag_error = 0;
            current_pos = zeros(1, num_round);
            current_path = survival_paths{current_root, num_round};
            current_remaining_roots_elimination(current_remaining_roots_elimination==current_root) = [];
            for current_round = 1 : num_round                
                remaining_message_current = remaining_message{current_round};
                corresp_pos = find(remaining_message_current==current_path(current_round));
                size_pos = length(corresp_pos);
                if (size_pos==0)
                    flag_error = 1;
                    break;
                end
                current_pos(current_round) = corresp_pos(1);
            end
            if (~flag_error) % if success, record it and remove it from the remaining set
                num_decoded_users = num_decoded_users + 1;
                final_decoded_users{num_decoded_users} = current_path;
                for current_round = 1 : num_round
                    remaining_message_current = remaining_message{current_round};
                    remaining_message_current(current_pos(current_round)) = [];
                    remaining_message{current_round} = remaining_message_current;
                end
            end
            
          else
              num_remainings = num_remainings + 1;
              remaining_roots(num_remainings) = current_root; % need to be further investigated, more than one path
          end
        end
        
        
        if (num_remainings==num_remainings_last)
          break;
        end

        num_remainings_last = num_remainings;
        remaining_roots_elimination = current_remaining_roots_elimination;

        % peel-off procedure
        for current_remain_root = remaining_roots
            current_possible_message = survival_paths{current_remain_root, num_round};
            for current_round = 2 : num_round
                current_remaining_messages = remaining_message{current_round};
                current_column = current_possible_message(:, current_round);
                if (isempty(current_column))
                    break;
                end
                num_candidate_current_round = length(current_column);
                num_elimination = 0;
                eliminated_messages = [];
                for l = 1 : num_candidate_current_round
                    if (sum(current_remaining_messages == current_column(l))==0)
                        %                pos_ = find(current_remaining_messages == current_column(l));
                        num_elimination = num_elimination + 1;
                        eliminated_messages(num_elimination) = l;
                        %                current_possible_message(l, :) = [];
                    end
                end
                current_possible_message(eliminated_messages, :) = [];
            end
            
            [num_survival_paths(current_remain_root, num_round), ~] = size(current_possible_message);
            survival_paths{current_remain_root, num_round} = current_possible_message;

        end
        
   end
   
   flag_resolve_repeated_root = 1;
   % repeated-root case, i.e, N trees with the same root have N survival paths.
   num_duplicate_roots = length(duplicate_decoded_roots);
   for num_group = 1 : num_duplicate_roots
       current_roots = duplicate_decoded_roots{num_group};
       size_current_roots = length(current_roots);
       current_survival_paths = survival_paths{current_roots(1), num_round};
       if (isempty(current_survival_paths))
           continue;
       end
       second_blocks = current_survival_paths(:, 2);
       size_second_blocks = length(second_blocks);
       if (size_second_blocks < 2)
           continue;
       end
       
       [~, unique_second_blocks_idx] = unique(second_blocks, 'stable');
       repeated_sets = [];
       repeated_idx = setdiff(1:size_second_blocks, unique_second_blocks_idx);
       for r_idx = repeated_idx
           delta_set = find(second_blocks==second_blocks(r_idx));
           repeated_sets = [repeated_sets, delta_set'];
       end
       if (~isempty(repeated_sets))
           repeated_sets = unique(repeated_sets, 'stable');
       end
       appear_once_sets = setdiff(1:size_second_blocks, repeated_sets);
       
       if (~isempty(appear_once_sets))
           size_once_appear = length(appear_once_sets);
           flag_resolve_repeated_root = 0;
           idx_r = 0;
           
           if (size_once_appear > size_current_roots)
               size_once_appear = size_current_roots; % just pick up the first size_current_roots paths.
           end
           
           for root_current = current_roots(1:size_once_appear)
               idx_r = idx_r + 1;
               num_survival_paths(root_current, num_round) = 1;
               survival_paths{root_current, num_round} = current_survival_paths(appear_once_sets(idx_r), :);
           end           
           for root_current = current_roots(size_once_appear+1:size_current_roots)
               num_survival_paths(root_current, num_round) = size_second_blocks-size_once_appear;
               survival_paths{root_current, num_round} = current_survival_paths(repeated_sets, :);
           end           
       end       
   end
   
   if (flag_resolve_repeated_root) % didn't resolve any repeated-root case;
       break;
   end
   
end

num_remains(iteration) = num_remainings;
end

PUPE = sum(num_remains) / num_iterations;
average_num_checks = sum(average_checks) / num_iterations;
end

注：

1、Markov superposition transmission (马尔科夫叠加传输):

2.2 Alex_Fengler_tree_codes

2.2.1 outerEncoder.m

function [ sparc_messages, C ] = outerEncoder( message_matrix, B, L, data_lengths, C )
%OUTERENCODER Summary of this function goes here
%   Detailed explanation goes here
    parity_lengths = log2(B) - data_lengths;% 计算每个符号的冗余位长度
    if nargin < 5 % 如果输入参数的个数小于 5（即没有提供 C），则调用 createParityMatrices 函数
        C = createParityMatrices(log2(B),parity_lengths);% 创建一个新的冗余校验矩阵 C (createParityMatrices定义见下)
    end
    sparc_messages  = encode(message_matrix, C, L, parity_lengths);% 调用 encode 函数，对输入的 message_matrix 进行编码。这个函数接受信息矩阵、冗余校验矩阵 C、轮数 L 以及冗余位长度，并返回编码后的消息矩阵 (encode函数定义见下)  
    sparc_messages = sparc_messages + 1;% 将编码后的消息矩阵的所有元素加 1(不清楚为什么加1)
end

function parity_m = createParityMatrices(b,parity_lengths)
    L = length(parity_lengths);
    data_lengths = b - parity_lengths;% 计算信息位的位数
    for i = 1:L %遍历生成
        data_bits   = i*b - sum(parity_lengths(1:i)); 
        parity_bits = parity_lengths(i);
        for k = 1:i
            parity_m{i}(:,k) = randi([0,2^data_lengths(k) - 1],parity_bits, 1);% 生成一个大小为 (parity_bits, 1) 的随机列向量，该列向量包含 [0, 2^data_lengths(k) - 1] 范围内的随机整数。这个向量用来存储奇偶校验位
        end
    end
end



function coded_matrix = encode(message_matrix, C, L, parity_lengths)
    K_a             = size(message_matrix,2);% 消息矩阵的列数
    coded_matrix    = zeros(L, K_a);% L为轮数
    
    for i = 1:L
        data_symbols    = message_matrix(1:i,:);% 包含消息矩阵的前 i 行
        parity_symbols  = zeros(1,K_a);% 用来表示当前奇偶校验比特的编码结果
        for k = 1:parity_lengths(i)% 遍历每个消息块，计算每个消息块对应的奇偶校验比特
            for j = 1:K_a 
                parity_bit   = LUT(bitand(C{i}(k,:)',data_symbols(:,j)));% 通过位与操作 bitand 计算奇偶校验比特的值。C{i}(k,:)' 表示当前奇偶校验比特的生成矩阵的第 k 行，data_symbols(:,j) 表示消息矩阵的第 j 列
                parity_symbols(j)  = parity_symbols(j) + parity_bit.*2^(parity_lengths(i));% 将计算得到的奇偶校验比特乘以2^(parity_lengths(i)) 实现左移操作，然后加到 parity_symbols 的相应位置。
                parity_symbols(j)  = bitshift(parity_symbols(j),-1);% 对 parity_symbols 中的当前元素进行右移操作，以便为下一轮计算腾出位置
            end

            
        end
        
        coded_matrix(i,:) = bitshift(data_symbols(i,:), parity_lengths(i)) + parity_symbols;% 通过左移和加法操作将 data_symbols 和 parity_symbols 合并得到当前轮的编码结果，即拼接信息位和校验位，并存储在 coded_matrix 中
    end
end

% Lookuptable for sums of bits
% Input:  list of integer numbers
% Output: 0 if the sum of the bit represented is even
%         1 if its odd
function out = LUT(symbol_array)
    out = mod(sum(sum(dec2bin(symbol_array) - '0')),2);
end% 功能同上LUT2函数，只是没有查表的操作


% 将矩阵中的二进制表示转换为整数
function out = bit2dec(bit_matrix)
    b   = size(bit_matrix,1);
    K   = size(bit_matrix,2);
    out = zeros(1,K);
    for i = 1:b
        out = out + bit_matrix(i,:).*2^(i-1);% 将当前行的二进制表示转换为整数，并加到结果向量 out 中
    end
end

注：

1、nargin函数：函数输入参数个数

2、randi函数：生成服从随机分布的随机数

2.2.2 outerDecoder.m

function [list, num_checks] = outerDecoder(symbol_list, C, parity_lengths, B, LT)
    % 用于从符号列表和相应的奇偶校验矩阵中生成消息列表 
    candidate_paths = symbol_list{1};% candidate_paths 的每列代表一个路径，每行代表路径上的一个节点，包括已经传输的符号。初始时，candidate_paths 包含的是第一轮传输的符号。
    L               = length(parity_lengths);% 奇偶校验码的长度
    num_checks = 0;% 检查的数量
    for l = 2:L% 循环遍历奇偶校验长度
        symbols = symbol_list{l};
        new_candidates  = [];
        for i_path = 1:size(candidate_paths,2)
            path = candidate_paths(1:l-1,i_path);% 遍历候选路径，取出当前路径
            for i_symbol = 1:length(symbols)
                symbol          = symbols(i_symbol);
                parity_symbol   = rem(symbol,2^parity_lengths(l));
                data_symbol     = bitshift(symbol, -parity_lengths(l));
                p               = parity_check_last([path;data_symbol], parity_symbol,C,log2(B),parity_lengths, LT);% 对于当前路径和当前符号，将符号分解成奇偶校验部分和数据部分，进行奇偶校验检查
                % parity_check_last为函数，定义见下
                num_checks = num_checks + 1;
                if p ==0 % 奇偶校验通过
                    new_candidates = [new_candidates [path;data_symbol]];% 将当前路径与符号的数据部分组合，并将其添加到新的候选路径列表中
                    if size(new_candidates,2) > 4000 % 检查新的候选路径列表的大小是否超过 4000
                        list = [];% 超过4000就清空结果列表
                        disp('Too many paths..');% 命令行中显示警告消息，表示路径过多
                        return;
                    end
                end
            end % 结束对当前符号的处理 
        end% 结束对当前路径的处理
        candidate_paths = new_candidates;% 更新 candidate_paths
        n_paths = size(candidate_paths,2);% 获取当前候选路径列表的大小
        
    end
    list = candidate_paths;
end


% 仅检查最后一个奇偶校验部分
function out = parity_check_last(data_symbols, parity_symbol, C, ~, parity_lengths, LT)
    L               = length(data_symbols);
    parity_lengths  = parity_lengths(1:L);% parity_lengths 包含了每个奇偶校验矩阵的奇偶校验位数量
    out             = 0;
    for i = 1:parity_lengths(L)% 循环遍历当前奇偶校验矩阵的奇偶校验位
        parity_bit = LUT2(bitand(C{L}(i,:)',data_symbols), LT);% 计算奇偶校验位，通过对数据符号和奇偶校验矩阵中的位进行按位与操作
        if rem(parity_symbol,2)~= parity_bit% 奇偶校验检查不通过
            out = 1;
            break;% 将 out 设置为 1，并跳出循环
        end
        parity_symbol = bitshift(parity_symbol,-1);% 将奇偶校验符号右移一位，准备进行下一位的比较
    end

end

% 用于位求和的查找表函数
% 输入参数是整数数组
% 输出值为 0，如果位求和为偶数
%         1 位求和为奇数
function out = LUT2(symbol_array, LT)
    binary = LT(symbol_array+1);% 使用输入数组 symbol_array 中的元素作为下标，从查找表 LT 中获取相应的二进制表示(matlab中下标从1开始)
    out = mod(sum(sum(binary)),2);% 第一次sum是对列求和，得出一个行向量，第二个sum是对行求和，得到最终的标量，然后模2判断奇偶
end

注：

1、size函数:返回数组大小

2、rem函数：返回除法的余数

3、bitshift函数：位左移或右移

4、bitand函数:按位与

2.2.3 UMAC_lossless_stitching_tree_codes.m

function [PUPE, average_num_checks] = UMAC_lossless_stitching_tree_codes(num_active, num_round, M, r,  num_iterations)
%  function [p_md,p_fa] = unsourcedSPARC(L, S, J, K_a, M, data_profile, P, iter, fading)
%unsourcedSPARC Compute the error rate of the MU-SPARC coding scheme
%
%   num_round:              Number of slots
%   num_active:            active users
%   M:              num of bits per section
%   r:   vector of the number of parity bits per section
%   num_iterations:           number of simulation runs
%
%   struct fading with the following fields
%   'type' = {'no_fading'(default),'uniform','exp','pathloss','shadowing_pathloss'}
%   'lower_limit'(default = 10),'upper_limit'(default = 30): limits for the uniform case

J = log2(M);
N               = 2^J*num_round;
data_profile  = J - r;


% Lookup table for outer code
LT  = createLT(2^max(data_profile));

num_checks = zeros(1, num_iterations);
p_md = 0;
p_fa = 0;

for iter=1:num_iterations
    
    % random messages
    msg = generateMessages(data_profile, num_active);

    % C is the matrix of parity checks needed for the decoder
    [symbol, C] = outerEncoder( msg, 2^J, num_round, data_profile );

    for num_r = 1 : num_round
        symbol_list{num_r} = symbol(num_r, :)-1;
    end
    
    [msg_list, num_checks(iter)]    = outerDecoder(symbol_list, C, r, 2^J, LT);

    [fn,fp]     = countErrors(msg_list, msg);
    p_md        = p_md + fn/num_active;
    p_fa        = p_fa + fp/length(msg_list);
end
p_md = p_md/num_iterations;
p_fa = p_fa/num_iterations;
PUPE = p_md + p_fa;
average_num_checks = sum(num_checks) / num_iterations;

end



%create lookuptable for sums of bits
function LT = createLT(N)
    LT = zeros(N,1);
    for i = 1:N
        LT(i) = mod(sum(dec2bin(i-1)-'0'),2);
    end
end

% Return: LxK matrix of symbols where each symbol has data_profile(l) bits
function msg_matrix = generateMessages(data_profile, K)
    L           = length(data_profile);
    msg_matrix  = zeros(L, K);
    for i = 1:L
        msg_matrix(i,:) = randi([0,2^data_profile(i) - 1], 1, K);
    end
end


% Create list of position indices from a support vector
function symbols = getSymbolList(rpos_matrix)
    L = size(rpos_matrix,2);
    for i = 1:L
        % outer decoder wants symbols from 0:B-1;
        symbols{i} = find(rpos_matrix(:,i)==1)' - 1;
    end
end

% Count the messages which miss in the output list, and the messages
% which were not transmitted but appeared in the output list
function [errs, false_positives] = countErrors(est, sparc_matrix)
    K       = size(est,2);
    errs    = size(sparc_matrix,2);
    
    false_positives = 0;
    for i = 1:K
        codeword = est(:,i);
        [~,indx]=ismember(codeword',sparc_matrix','rows');
        if (indx~=0)
            sparc_matrix(:,indx) = [];
            errs = errs -1;
        else
            false_positives = false_positives + 1;
        end        
    end
end

三、CRC-aided SPARCs for URA

3.1 CRC_check.m

function check = CRC_check(r_c, g)

check = 0;
[~, r] = gfdeconv(fliplr(r_c), fliplr(g));
if (sum(r)==0)
    check = 1;
end

end
%同2.1.1

3.2 outer_CRC_checks_BMST_efficient.m

function [current_blocks, check] = outer_CRC_checks_BMST_efficient(previous_blocks, previous_superimposed_sum, current_path, current_block, num_protect_section, m, memory, r, CRC_poly)
% the (efficient) CRC check procedure starting with substracting BMST;
% store previous original_block_bits so that the same info doesn't need to compute twice.
% start from the second layer.

% current_blocks: the recovered info bits up to current blocks;
% check: indicator regarding whether the current bits satisfy the CRC condtion;

% previous_blocks: all previous recovered info bits;
% previous_superimposed_sum: the previoused computed bit-wise sum;
% current_path: the current path including all nodes;
% current_block: the current stage or round;
% num_protect_section: the number of protected sections;
% m: the number of bits for each block, i.e., log2(M);
% memory: the memory parameter used in BMST encoding;
% r: a vector indicating the number of redundant bits at each round;
% CRC_poly: the chosen CRC generator polynomial;

    check = 0;
    original_block_bits = zeros(current_block, m);
    original_block_bits(1:current_block-1, :) = previous_blocks;


    current_block_bits = de2bi(current_path(current_block)-1, m);
    original_block_bits(current_block, :) = mod(current_block_bits-previous_superimposed_sum, 2);

    current_path_info_bits = [];
    start_cc = 0;

    for cc = current_block - num_protect_section+1 : current_block-1
      end_cc = start_cc + m - r(cc);
      current_path_info_bits(start_cc+1: end_cc) = original_block_bits(cc, 1:m-r(cc));
      start_cc = end_cc;
    end


    current_blocks = original_block_bits;

    check_bits = [current_path_info_bits, original_block_bits(current_block, :)];
    check = CRC_check(check_bits, CRC_poly);

end

%同2.1.2

3.3 CRC_encoding.m

function c = CRC_encoding(m, g)

% m : information vector
% g : generator polynomial
% c : the CRC part
% since the order in gfdeconv is ascending, we need do an "fliplr"
% operation.

c = zeros(1, length(g)-1);
m_r = [m zeros(1, length(g)-1)];
[~, r] = gfdeconv(fliplr(m_r), fliplr(g));
c(1:length(r)) = r;
c = fliplr(c);

end
%同2.1.3

3.4 outer_encoding_BMST.m

function outer_encoded = outer_encoding_BMST(message, memory, m, num_round, r, protect_sections, num_info_bits_round, orderings)
  % the CRC-BMST encoding procedure
  % we only protect information bits via CRC codes

  % message: message bits;
  % memory: the encoding memory parameter used in BMST encoding;
  % m: the number of bits for each section;
  % num_round: the original messages transmitted by each user are splitted into such rounds;
  % r: a vector indicating the number of redundant bits at each round;
  % protect_sections: the number of preceding sections protected at current round;
  % num_info_bits_round: a vector indicating the number of info bits at each round;
  % orderings: interleavers used in CRC-BMST encoding procedure;


% the CRC generator polynomials with different number of CRC bits
poly{4} = [1 0 0 1 1]; % the number of info bits up to length 11

%poly{5} = [1 0 0 1 0 1]; % up to length 26
poly{5} = [1 0 1 0 1 1]; % up to length 10

%poly{6} = [1 0 0 0 0 1 1]; % up to length 57
poly{6} = [1 0 0 0 1 1 1]; % up to (information) length 25

poly{7} = [1 0 1 1 0 1 1 1]; % up to length 56

%poly{8} = [1 0 0 1 0 1 1 1 1]; % up to length 119
%poly{8} = [1 0 0 0 0 0 1 1 1]; % up to length 119
poly{8} = [1 1 1 0 1 0 1 1 1]; % up to length 9

%poly{9} = [1 0 1 0 0 1 0 1 1 1]; % up to length 246
%poly{9} = [1 0 1 1 1 1 1 0 1 1]; % up to length 246
%poly{9} = [1 1 0 0 0 0 1 0 1 1]; % up to length 13
poly{9} = [1 0 0 1 1 1 1 0 0 1]; % up to length 8

%poly{10} = [1 1 0 0 0 1 1 0 0 1 1]; % up to length 501
%poly{10} = [1 0 1 0 1 1 1 0 0 1 1]; % up to length 21
%poly{10} = [1 0 1 0 0 0 1 1 1 0 1]; % up to length 12
poly{10} = [1 0 1 0 0 1 1 0 1 1 1]; % up to length 5

% poly{11} = [1 0 0 1 1 1 1 0 1 0 1 1]; % up to length 4
poly{11} = [1 0 1 0 1 1 1 0 0 0 1 1]; % up to length 12

poly{12} = [1 0 1 0 0 1 0 0 1 1 1 1 1]; % up to length 11

poly{13} = [1 0 0 0 0 1 0 1 1 0 1 1 1 1]; % up to length 11

poly{14} = [1 0 0 0 1 1 0 1 1 1 0 0 0 1 1]; % up to length 11


  outer_encoded(1:m) = message(1:m);
  message_last_end = m;
  original_encoded_block = zeros(num_round, m);
  original_encoded_block(1, :) = message(1:m);

  for current_round = 2 : num_round
    protected_info_bits = [];
    num_protect_section = protect_sections(current_round);
    encoded_last_end = (current_round-1)*m;
    info_current = message(message_last_end+1 : message_last_end + m -r(current_round));
    start_ = 0;
    if (current_round > num_protect_section)
      start_ = num_info_bits_round(current_round-num_protect_section);
    end
    protected_info_bits =message(start_+1 : num_info_bits_round(current_round));
    outer_encoded_current = [protected_info_bits, CRC_encoding(protected_info_bits, poly{r(current_round)})];


    len_outer_current = length(outer_encoded_current);
    start_block_BMST = max(current_round-memory, 1);
    outer_encoded_block = outer_encoded_current(len_outer_current-m+1:len_outer_current);
    original_encoded_block(current_round, :) = outer_encoded_block;
    num_m = 0;
    for sb = current_round-1 : -1 : start_block_BMST
      num_m = num_m + 1;
      permuting_block = original_encoded_block(sb, :);
      outer_encoded_block = mod(outer_encoded_block + permuting_block(orderings(num_m, :)),2);
    end
    outer_encoded(encoded_last_end+1:encoded_last_end+m) = outer_encoded_block;
    message_last_end = message_last_end + m - r(current_round);
  end
  
end
%同2.1.4

3.5 Outer_UMAC_stitching_CRC_BMST.m

function [final_decoded_users] = Outer_UMAC_stitching_CRC_BMST(num_candidates,orderings, outer_decoded_pos, num_round, M, r, protect_sections, memory, poly, duplicate_decoded_roots)
  % tree decoder for UMAC, i.e., stitching procedure
  % only output unordered list of decoded messages transmitted over unsourced MAC.

  % final_decoded_users: the unordered list of decoded messages;

  % num_active: the number of active users;
  % num_round:  the number of chunks;
  % M: the size of codebook;
  % n: the length of codewords;
  % r: a vector consisted with entries showing the number of redundant bits for each chunk;
  % protect_sections: indicate the CRC codeword protecting how many previously consecutive chunks;
  % memory: the memory parameter used in BMST-CRC encoding;
  % poly: CRC generator polynomial;
  % duplicate_decoded_roots: duplicated decoded roots for different decoding trees;

  
  m = log2(M);
  final_decoded_users = [];

   % stitching procedure

   decoded_CRC_encoded_pos = outer_decoded_pos;

   % with the aid of CRC codes, record the survival paths
   
   num_survival_paths = zeros(num_candidates, num_round);
   survival_current_blocks = cell(num_candidates, num_candidates);

   survival_paths = cell(num_candidates, num_round);  % each cell represents a matrix containing all survival paths.
   
   for current_user = 1 : num_candidates
     survival_paths{current_user, 1} = decoded_CRC_encoded_pos(1, current_user);
     num_survival_paths(current_user, 1) = 1;
     survival_current_blocks{current_user, 1} = de2bi(decoded_CRC_encoded_pos(1, current_user)-1, m);
   end


   for current_round = 2 : num_round
       survival_blocks_prev = cell(num_candidates, num_candidates);
       survival_blocks_prev = survival_current_blocks;
       survival_current_blocks = cell(num_candidates, num_candidates);
       for current_root = 1 : num_candidates
           num_paths_current = 0;
           num_paths_prev = num_survival_paths(current_root, current_round-1);
           survival_paths_prev = survival_paths{current_root, current_round-1};
           survival_path_current = [];
           unique_msg_current = unique(decoded_CRC_encoded_pos(current_round, :), "stable");
           num_protect_section = protect_sections(current_round);
           for current_path_num = 1 : num_paths_prev
               survival_path_prev = survival_paths_prev(current_path_num, :);
               previous_blocks = survival_blocks_prev{current_root, current_path_num};
               superimposed_sum = zeros(1, m);
               num_m = 0;
               start_block_BMST_ = max(current_round-memory, 1);
               for prev_batch = current_round-1 : -1 : start_block_BMST_
                   num_m = num_m + 1;
                   superimposed_sum = superimposed_sum + previous_blocks(prev_batch, orderings(num_m,:));
               end
               
               previous_superimposed_sum = mod(superimposed_sum, 2);
               for current_pos_considered = unique_msg_current
                   CRC_poly = poly{r(current_round)};
                   current_path = [survival_path_prev, current_pos_considered];
                   current_blocks = [];
                   [current_blocks,check] = outer_CRC_checks_BMST_efficient(previous_blocks, previous_superimposed_sum, current_path, current_round, num_protect_section, m, memory, r, CRC_poly);
                   if (check)
                       num_paths_current = num_paths_current + 1;
                       survival_path_current(num_paths_current, 1:current_round) = [survival_path_prev, current_pos_considered];
                       survival_current_blocks{current_root, num_paths_current} = current_blocks;
                   end
               end
           end
           num_survival_paths(current_root, current_round) = num_paths_current;
           survival_paths{current_root, current_round} = survival_path_current;
       end
       %         disp(current_round);
   end
   
   
   % peeling decoding procedure in the sense that peel off nodes included in the unique path
   remaining_message = cell(1, num_round);
   for current_round = 1 : num_round
     remaining_message{current_round} = decoded_CRC_encoded_pos(current_round, :);
   end

   num_remainings_last = num_candidates;
   remaining_roots_elimination = 1 : num_candidates;
   final_num_survival_paths = num_survival_paths(:, num_round);
   num_detected_error_trees = sum(final_num_survival_paths==0);
   remaining_roots_elimination(final_num_survival_paths==0) = []; % excluding the detected failure trees.
   num_decoded_users = 0;

while (1)     
        
   while (1)
        num_remainings = 0;
        remaining_roots = [];
        current_remaining_roots_elimination = remaining_roots_elimination;

        % eliminating the corrected messages
        for current_root  = remaining_roots_elimination
          if (num_survival_paths(current_root, num_round)==1)
            flag_error = 0;
            current_pos = zeros(1, num_round);
            current_path = survival_paths{current_root, num_round};
            current_remaining_roots_elimination(current_remaining_roots_elimination==current_root) = [];
            for current_round = 1 : num_round                
                remaining_message_current = remaining_message{current_round};
                corresp_pos = find(remaining_message_current==current_path(current_round));
                size_pos = length(corresp_pos);
                if (size_pos==0)
                    flag_error = 1;
                    break;

                end
                current_pos(current_round) = corresp_pos(1);
            end
            if (~flag_error) % if success, record it and remove it from the remaining set
                num_decoded_users = num_decoded_users + 1;
                final_decoded_users{num_decoded_users} = current_path;
                for current_round = 1 : num_round
                    remaining_message_current = remaining_message{current_round};
                    remaining_message_current(current_pos(current_round)) = [];
                    remaining_message{current_round} = remaining_message_current;
                end
            end
            
          else
              num_remainings = num_remainings + 1;
              remaining_roots(num_remainings) = current_root; % need to be further investigated, more than one path
          end
        end
        
        
        if (num_remainings==num_remainings_last)
          break;
        end

        num_remainings_last = num_remainings;
        remaining_roots_elimination = current_remaining_roots_elimination;

        % peel-off procedure
        for current_remain_root = remaining_roots
            current_possible_message = survival_paths{current_remain_root, num_round};
            for current_round = 2 : num_round
                current_remaining_messages = remaining_message{current_round};
                current_column = current_possible_message(:, current_round);
                if (isempty(current_column))
                    break;
                end
                num_candidate_current_round = length(current_column);
                num_elimination = 0;
                eliminated_messages = [];
                for l = 1 : num_candidate_current_round
                    if (sum(current_remaining_messages == current_column(l))==0)
                        num_elimination = num_elimination + 1;
                        eliminated_messages(num_elimination) = l;
                    end
                end
                current_possible_message(eliminated_messages, :) = [];
            end
            
            [num_survival_paths(current_remain_root, num_round), ~] = size(current_possible_message);
            survival_paths{current_remain_root, num_round} = current_possible_message;

        end       
   end
   
   flag_resolve_repeated_root = 1;

   % repeated-root case, i.e, N trees with the same root have N survival paths.
   num_duplicate_roots = length(duplicate_decoded_roots);
   for num_group = 1 : num_duplicate_roots
       current_roots = duplicate_decoded_roots{num_group};
       size_current_roots = length(current_roots);
       current_survival_paths = survival_paths{current_roots(1), num_round};
       if (isempty(current_survival_paths))
           continue;
       end
       second_blocks = current_survival_paths(:, 2);
       size_second_blocks = length(second_blocks);
       if (size_second_blocks < 2)
           continue;
       end
       
       [~, unique_second_blocks_idx] = unique(second_blocks, 'stable');
       repeated_sets = [];
       repeated_idx = setdiff(1:size_second_blocks, unique_second_blocks_idx);
       for r_idx = repeated_idx
           delta_set = find(second_blocks==second_blocks(r_idx));
           repeated_sets = [repeated_sets, delta_set'];
       end
       if (~isempty(repeated_sets))
           repeated_sets = unique(repeated_sets, 'stable');
       end
       appear_once_sets = setdiff(1:size_second_blocks, repeated_sets);
       
       if (~isempty(appear_once_sets))
           size_once_appear = length(appear_once_sets);
           flag_resolve_repeated_root = 0;
           idx_r = 0;
           if (size_once_appear > size_current_roots)
               size_once_appear = size_current_roots; % just pick up the first size_current_roots paths.
           end
           
           for root_current = current_roots(1:size_once_appear)
               idx_r = idx_r + 1;
               num_survival_paths(root_current, num_round) = 1;
               survival_paths{root_current, num_round} = current_survival_paths(appear_once_sets(idx_r), :);
           end           
           for root_current = current_roots(size_once_appear+1:size_current_roots)
               num_survival_paths(root_current, num_round) = size_second_blocks-size_once_appear;
               survival_paths{root_current, num_round} = current_survival_paths(repeated_sets, :);
           end           
       end       
   end
   
   if (flag_resolve_repeated_root) % didn't resolve any repeated-root case;
       break;
   end
   
end 

end

3.6 MAP_AMP_hybrid_UMAC_CRC_BMST.m

function [PUPE_MD, PUPE_FA] = MAP_AMP_hybrid_UMAC_CRC_BMST(num_active, num_round, num_info_bits, M, n, P, r, protect_sections, num_trials, memory)
  % encoding and decoding of the proposed concatenated coding system for
  % BMST-CRC codes and SPARCs.
  % only output unordered list of decoded messages transmitted over unsourced MAC.
  
  % PUPE_MD: PUPE for miss detection;
  % PUPE_FA: PUPE for false alarm;

  % num_active: the number of active users;
  % num_round:  the number of chunks;
  % num_info_bits: a vector consisted with entries showing the number of  info bits for each chunk;
  % M: the size of codebook;
  % n: the length of codewords;
  % P: the averag power constraint;
  % r: a vector consisted with entries showing the number of redundant bits for each chunk;
  % protect_sections: indicate the CRC codeword protecting how many previously consecutive chunks;
  % num_trials: the number of simulation runs;
  % memory: the memory parameter used in BMST-CRC encoding;

  
% the CRC generator polynomials with different number of CRC bits
  poly{4} = [1 0 0 1 1]; % number of info bits up to length 11

  % poly{5} = [1 0 0 1 0 1]; % up to length 26
  poly{5} = [1 0 1 0 1 1]; % up to length 10

  % poly{6} = [1 0 0 0 0 1 1]; % up to length 57
  poly{6} = [1 0 0 0 1 1 1]; % up to (information) length 25

  poly{7} = [1 0 1 1 0 1 1 1]; % up to length 56

  % poly{8} = [1 0 0 1 0 1 1 1 1]; % up to length 119
  % poly{8} = [1 0 0 0 0 0 1 1 1]; % up to length 119
  poly{8} = [1 1 1 0 1 0 1 1 1]; % up to length 9

  % poly{9} = [1 0 1 0 0 1 0 1 1 1]; % up to length 246
  % poly{9} = [1 0 1 1 1 1 1 0 1 1]; % up to length 246
  % poly{9} = [1 1 0 0 0 0 1 0 1 1]; % up to length 13
  poly{9} = [1 0 0 1 1 1 1 0 0 1]; % up to length 8

  % poly{10} = [1 1 0 0 0 1 1 0 0 1 1]; % up to length 501
  % poly{10} = [1 0 1 0 1 1 1 0 0 1 1]; % up to length 21
  % poly{10} = [1 0 1 0 0 0 1 1 1 0 1]; % up to length 12
  % poly{10} = [1 0 1 0 0 1 1 0 1 1 1]; % up to length 5

  poly{11} = [1 0 0 1 1 1 1 0 1 0 1 1]; % up to length 4 

L = num_active;
K = num_active;
m = log2(M);

P_l = P;

% randomly generate interleavers
orderings = zeros(memory, m);
inverse_orderings = zeros(memory, m);
for mmr = 1 : memory
  perm_current = randperm(m);
  inverse_perm_current(perm_current) = 1 : m;
  orderings(mmr, :) = perm_current;
  inverse_orderings(mmr, :) = inverse_perm_current;
end


% generate the measurement matrix A by randomly choosing n rows from Hadamard matrix.
n_2power = M;
ordering =  randperm(n_2power-1, n)+1; 


num_info_bits_round(1) = m;
for current_round = 2 : num_round
  num_info_bits_round(current_round) = num_info_bits_round(current_round-1) + m - r(current_round);
end

extra_candidates = 10;
num_candidates = K + extra_candidates;
num_missed_detection = zeros(num_trials, 1);
num_false_alarm = zeros(num_trials, 1);

for trial = 1 : num_trials
   rng('shuffle');
   outer_encoded = zeros(L, m*num_round);
   message = zeros(L, num_info_bits);

   % outer encoding via CRC codes at each round
   for l = 1 : L
     message(l, :) = randi(2, 1, num_info_bits)-1;
     if (sum(message(l, 1:m))==0) % to avoid taking the first all-one codeword
       message(l, m) = 1;
     end
     outer_encoded(l, :) = outer_encoding_BMST(message(l, :), memory, m, num_round, r, protect_sections, num_info_bits_round, orderings);
   end

   % inner encoding and decoding via SPARCs
   CRC_encoded_pos = zeros(num_round, L);

   decoded_CRC_encoded_candidates = zeros(num_round, num_candidates);

   for current_round = 1 : num_round
     beta = zeros(M, 1);
     x = zeros(n, 1);
     info_bits = zeros(m, L);
     info_pos = zeros(1, L);
     decoded_bits = zeros(m, L);
     last_end = (current_round-1)*m;
     
      % inner encoding via SPARCs chunk by chunk
     for l = 1 : L
       current_CRC_encoded_bits = outer_encoded(l, last_end+1:last_end+m);
       CRC_encoded_pos(current_round, l) = bi2de(current_CRC_encoded_bits)+1;
       beta(CRC_encoded_pos(current_round, l)) = beta(CRC_encoded_pos(current_round, l)) + 1;
     end
     
     x = sqrt(P)*FWHT_user_subsampling(1, beta, n_2power, M, ordering);
   
     z = randn(n, 1);  % additive Gaussian noise

     y = x + z; % the received signal

     % hybrid inner decoding chunk by chunk
     [decoded_CRC_encoded_candidates(current_round, :)] = MAP_AMP_Hybrid_SPARCs_UMAC(y, n, M, K, extra_candidates, ordering, P_l, 100);

   end
   
   
   
   % count duplicate roots
   decoded_roots_msg= decoded_CRC_encoded_candidates(1, :);

   num_decoded_duplicate_roots = 0;
   duplicate_decoded_roots = {};
   [~, unique_roots_idx] = unique(decoded_roots_msg, 'stable');
   roots_duplicate_idx = setdiff(1:num_candidates, unique_roots_idx);
   
   while (~isempty(roots_duplicate_idx))
       num_decoded_duplicate_roots = num_decoded_duplicate_roots + 1;
       current_root_idx = roots_duplicate_idx(1);
       current_root_msg = decoded_roots_msg(current_root_idx);
       duplicate_decoded_roots_current = find(decoded_roots_msg==current_root_msg);
       duplicate_decoded_roots{num_decoded_duplicate_roots} = duplicate_decoded_roots_current;
       roots_duplicate_idx = setdiff(roots_duplicate_idx, duplicate_decoded_roots_current);
   end

   
   % stitching procedure
   final_decoded_users = Outer_UMAC_stitching_CRC_BMST(num_candidates, orderings, decoded_CRC_encoded_candidates, num_round, M, r, protect_sections, memory, poly, duplicate_decoded_roots);
 
   num_decoded_users = size(final_decoded_users, 2);   
   decoded_roots = [];
   for c_u = 1 : num_decoded_users
       decoded_msg = final_decoded_users{c_u};
       decoded_roots(c_u) = decoded_msg(1);
   end

   user_roots = CRC_encoded_pos(1, :);
   num_correct_decoded = 0;
   
   % remove possible duplicates
   [~, unique_idx] = unique(decoded_roots, 'stable');
   root_duplicate_idx = setdiff(1:num_decoded_users, unique_idx);
   duplicate_set = [];
   current_user = [];
   for idx = 1 : length(root_duplicate_idx)
       current = root_duplicate_idx(idx);
       current_user = final_decoded_users{current};
       root_duplicates = setdiff(find(decoded_roots==decoded_roots(current)), current);
       for compare_roots = root_duplicates
           if (~ismember(compare_roots, duplicate_set))
               compare_user = final_decoded_users{compare_roots};
               if (prod(compare_user==current_user))
                   duplicate_set = [duplicate_set, compare_roots];
               end
           end
       end       
   end
   
   if (~isempty(duplicate_set))
       final_decoded_users(duplicate_set) = [];
       num_decoded_users = size(final_decoded_users, 2);
   end
   
   for c_u = 1 : num_decoded_users
       decoded_msg = final_decoded_users{c_u};
       decoded_user = find(user_roots == decoded_msg(1));
       
       if (decoded_user)
           if (sum(prod(CRC_encoded_pos(:, decoded_user)==decoded_msg'))) % paths with the same roots are included
               num_correct_decoded = num_correct_decoded + 1;
           end
       end       
   end
   
   num_missed_detection(trial) = num_active - num_correct_decoded;
   num_false_alarm(trial) = num_decoded_users - num_correct_decoded;

end

PUPE_MD = sum(num_missed_detection) / (num_trials*num_active);
PUPE_FA = sum(num_false_alarm) / (num_trials*num_active);

end

3.7 MAP_AMP_hybrid_UMAC_CRC_BMST_extra_iteration.m

function [PUPE_MD, PUPE_FA, PUPE_MD_again, PUPE_FA_again] = MAP_AMP_hybrid_UMAC_CRC_BMST_extra_iteration(num_active, num_round, num_info_bits, M, n, P, r, protect_sections, num_trials, memory)
  % encoding and decoding of the proposed concatenated coding system for
  % BMST-CRC codes and SPARCs with an extra iteration for fair comparison.

  % PUPE_MD: PUPE for miss detection;
  % PUPE_FA: PUPE for false alarm;
  % PUPE_MD_again: PUPE for miss detection after one extra iteration;
  % PUPE_FA_again: PUPE for false alarm after one extra iteration;

  % num_active: the number of active users;
  % num_round:  the number of chunks;
  % num_info_bits: a vector consisted with entries showing the number of  info bits for each chunk;
  % M: the size of codebook;
  % n: the length of codewords;
  % P: the averag power constraint;
  % r: a vector consisted with entries showing the number of redundant bits for each chunk;
  % protect_sections: indicate the CRC codeword protecting how many previously consecutive chunks;
  % num_trials: the number of simulation runs;
  % memory: the memory parameter used in BMST-CRC encoding;

% the CRC generator polynomials with different number of CRC bits
poly{4} = [1 0 0 1 1]; % up to length 11

% poly{5} = [1 0 0 1 0 1]; % up to length 26
poly{5} = [1 0 1 0 1 1]; % up to length 10

% poly{6} = [1 0 0 0 0 1 1]; % up to length 57
poly{6} = [1 0 0 0 1 1 1]; % up to (information) length 25

poly{7} = [1 0 1 1 0 1 1 1]; % up to length 56

% poly{8} = [1 0 0 1 0 1 1 1 1]; % up to length 119
% poly{8} = [1 0 0 0 0 0 1 1 1]; % up to length 119
poly{8} = [1 1 1 0 1 0 1 1 1]; % up to length 9

% poly{9} = [1 0 1 0 0 1 0 1 1 1]; % up to length 246
% poly{9} = [1 0 1 1 1 1 1 0 1 1]; % up to length 246
% poly{9} = [1 1 0 0 0 0 1 0 1 1]; % up to length 13
poly{9} = [1 0 0 1 1 1 1 0 0 1]; % up to length 8

% poly{10} = [1 1 0 0 0 1 1 0 0 1 1]; % up to length 501
% poly{10} = [1 0 1 0 1 1 1 0 0 1 1]; % up to length 21
% poly{10} = [1 0 1 0 0 0 1 1 1 0 1]; % up to length 12
% poly{10} = [1 0 1 0 0 1 1 0 1 1 1]; % up to length 5

% poly{11} = [1 0 0 1 1 1 1 0 1 0 1 1]; % up to length 4
poly{11} = [1 0 1 0 1 1 1 0 0 0 1 1]; % up to length 12

poly{12} = [1 0 1 0 0 1 0 0 1 1 1 1 1]; % up to length 11

poly{13} = [1 0 0 0 0 1 0 1 1 0 1 1 1 1]; % up to length 11

poly{14} = [1 0 0 0 1 1 0 1 1 1 0 0 0 1 1]; % up to length 11

L = num_active;
K = num_active;
m = log2(M);

P_l = P;

% randomly generate interleavers
orderings = zeros(memory, m);
inverse_orderings = zeros(memory, m);
for mmr = 1 : memory
  perm_current = randperm(m);
  inverse_perm_current(perm_current) = 1 : m;
  orderings(mmr, :) = perm_current;
  inverse_orderings(mmr, :) = inverse_perm_current;
end


% generate the measurement matrix A by randomly choosing n rows from Hadamard matrix.
n_2power = M;
ordering =  randperm(n_2power-1, n)+1; 

num_info_bits_round(1) = m;

for current_round = 2 : num_round
  num_info_bits_round(current_round) = num_info_bits_round(current_round-1) + m - r(current_round);
end

extra_candidates = 10;
num_candidates = K + extra_candidates;
num_missed_detection_again = zeros(num_trials, 1);
num_false_alarm_again = zeros(num_trials, 1);
num_missed_detection = zeros(num_trials, 1);
num_false_alarm = zeros(num_trials, 1);


for trial = 1 : num_trials
   rng('shuffle');
   outer_encoded = zeros(L, m*num_round);
   message = zeros(L, num_info_bits);

   % outer encoding via CRC codes at each round
   for l = 1 : L
     message(l, :) = randi(2, 1, num_info_bits)-1;
     if (sum(message(l, 1:m))==0) % to avoid taking the first all-one codeword
       message(l, m) = 1;
     end
     outer_encoded(l, :) = outer_encoding_BMST(message(l, :), memory, m, num_round, r, protect_sections, num_info_bits_round, orderings);
   end

   % inner encoding and decoding via SPARCs
   CRC_encoded_pos = zeros(num_round, L);

   decoded_CRC_encoded_candidates = zeros(num_round, num_candidates);

  
   y = zeros(n, num_round);
   
   for current_round = 1 : num_round
     beta = zeros(M, 1);
     last_end = (current_round-1)*m;
     
      % inner encoding via SPARCs patch by patch
     for l = 1 : L
       current_CRC_encoded_bits = outer_encoded(l, last_end+1:last_end+m);
       CRC_encoded_pos(current_round, l) = bi2de(current_CRC_encoded_bits)+1;
       beta(CRC_encoded_pos(current_round, l)) = beta(CRC_encoded_pos(current_round, l)) + 1;
     end
     
     
     x = sqrt(P)*FWHT_user_subsampling(1, beta, n_2power, M, ordering);
    
     z = randn(n, 1);  
     y(:, current_round) = x + z;

     % MAP+AMP hybrid (inner) decoding pitch by pitch
     [decoded_CRC_encoded_candidates(current_round, :)] = MAP_AMP_Hybrid_SPARCs_UMAC(y(:, current_round), n, M, K, extra_candidates, ordering, P_l, 100);

   end
   
      
   % count duplicate roots
   decoded_roots_msg= decoded_CRC_encoded_candidates(1, :);

   num_decoded_duplicate_roots = 0;
   duplicate_decoded_roots = {};
   [~, unique_roots_idx] = unique(decoded_roots_msg, 'stable');
   roots_duplicate_idx = setdiff(1:num_candidates, unique_roots_idx);
   
   while (~isempty(roots_duplicate_idx))
       num_decoded_duplicate_roots = num_decoded_duplicate_roots + 1;
       current_root_idx = roots_duplicate_idx(1);
       current_root_msg = decoded_roots_msg(current_root_idx);
       duplicate_decoded_roots_current = find(decoded_roots_msg==current_root_msg);
       duplicate_decoded_roots{num_decoded_duplicate_roots} = duplicate_decoded_roots_current;
       roots_duplicate_idx = setdiff(roots_duplicate_idx, duplicate_decoded_roots_current);
   end   
   
   % stitching procedure
   final_decoded_users = Outer_UMAC_stitching_CRC_BMST(num_candidates, orderings, decoded_CRC_encoded_candidates, num_round, M, r, protect_sections, memory, poly, duplicate_decoded_roots);

   
   num_decoded_users = size(final_decoded_users, 2);
   
   for c_u = 1 : num_decoded_users
       decoded_msg = final_decoded_users{c_u};
       decoded_roots(c_u) = decoded_msg(1);
   end

   
   % remove possible duplicates
   [~, unique_idx] = unique(decoded_roots, 'stable');
   root_duplicate_idx = setdiff(1:num_decoded_users, unique_idx);
   duplicate_set = [];
   for idx = 1 : length(root_duplicate_idx)
       current = root_duplicate_idx(idx);
       current_user = final_decoded_users{current};
       root_duplicates = setdiff(find(decoded_roots==decoded_roots(current)), current);
       for compare_roots = root_duplicates
           if (~ismember(compare_roots, duplicate_set))
               compare_user = final_decoded_users{compare_roots};
               if (prod(compare_user==current_user))
                   duplicate_set = [duplicate_set, compare_roots];
               end
           end
       end       
   end
   
   if (~isempty(duplicate_set))
       final_decoded_users(duplicate_set) = []; % remove the duplicated copy
       num_decoded_users = size(final_decoded_users, 2);
   end
   
   num_correct_decoded = 0;
   user_roots = CRC_encoded_pos(1, :);
   decoded_msgs = zeros(num_decoded_users, num_round);
   for c_u = 1 : num_decoded_users
       decoded_msg = final_decoded_users{c_u};
       decoded_msgs(c_u, :) = decoded_msg;
       decoded_user = find(user_roots == decoded_msg(1));
       if (decoded_user)
           if (sum(prod(CRC_encoded_pos(:, decoded_user)==decoded_msg'))) % paths with the same roots are included
               num_correct_decoded = num_correct_decoded + 1;
           end
       end
   end
   
   num_missed_detection(trial) = num_active - num_correct_decoded;
   num_false_alarm(trial) = num_decoded_users - num_correct_decoded;
   

   % one extra iteration, i.e., subtract the decoded users' messages and
   % run our proposed decoding algorithm for the remaining;

   num_left = K-num_decoded_users;
   if (num_left ==0)
       num_missed_detection_again(trial)  = num_missed_detection(trial);
       num_false_alarm_again(trial) = num_false_alarm(trial);
       continue;
   end
   
   num_left_candidates = num_left + extra_candidates;
   decoded_CRC_encoded_candidates_again = zeros(num_round, num_left_candidates);
   decoded_CRC_encoded_candidates_again_MRC = zeros(num_round, num_left_candidates);
   for current_round = 1 : num_round
       y_current = y(:, current_round);
       decoded_pos_current = decoded_msgs(:, current_round);
       beta = zeros(M, 1);
       for pos_idx = 1 : num_decoded_users
           beta(decoded_pos_current(pos_idx)) = beta(decoded_pos_current(pos_idx)) + 1;
       end
       x = sqrt(P)*FWHT_user_subsampling(1, beta, n_2power, M, ordering);
       y_current = y_current - x;
       [decoded_CRC_encoded_candidates_again(current_round, :)] = AMP_SPARCs_UMAC(y_current, n, M, num_left, extra_candidates, ordering, P_l, 100);
       [~, decoded_CRC_encoded_candidates_again_MRC(current_round, :)] = maxk(FWHT_user_subsampling(2, y_current, n_2power, M, ordering), num_left_candidates);
   end
   
   % count duplicate roots
   decoded_roots_msg= decoded_CRC_encoded_candidates_again(1, :);

   num_decoded_duplicate_roots = 0;
   duplicate_decoded_roots = {};
   [~, unique_roots_idx] = unique(decoded_roots_msg, 'stable');
   roots_duplicate_idx = setdiff(1:num_left_candidates, unique_roots_idx);
   
   while (~isempty(roots_duplicate_idx))
       num_decoded_duplicate_roots = num_decoded_duplicate_roots + 1;
       current_root_idx = roots_duplicate_idx(1);
       current_root_msg = decoded_roots_msg(current_root_idx);
       duplicate_decoded_roots_current = find(decoded_roots_msg==current_root_msg);
       duplicate_decoded_roots{num_decoded_duplicate_roots} = duplicate_decoded_roots_current;
       roots_duplicate_idx = setdiff(roots_duplicate_idx, duplicate_decoded_roots_current);
   end   
   
   
   % stitching procedure again
   final_decoded_users_again = Outer_UMAC_stitching_CRC_BMST(num_left_candidates, orderings, decoded_CRC_encoded_candidates_again, num_round, M, r, protect_sections, memory, poly, duplicate_decoded_roots);
     
   num_decoded_users_again = size(final_decoded_users_again, 2);
   
   if (num_decoded_users_again == 0)
       num_missed_detection_again(trial)  = num_missed_detection(trial);
       num_false_alarm_again(trial) = num_false_alarm(trial);
       continue;
   end 
   
   
   decoded_roots = [];
   for c_u = 1 : num_decoded_users_again
       decoded_msg = final_decoded_users_again{c_u};
       decoded_roots(c_u) = decoded_msg(1);
   end

   user_roots = CRC_encoded_pos(1, :);
   
   % remove possible duplicates
   [~, unique_idx] = unique(decoded_roots, 'stable');
   root_duplicate_idx = setdiff(1:num_decoded_users_again, unique_idx);
   duplicate_set = [];
   for idx = 1 : length(root_duplicate_idx)
       current = root_duplicate_idx(idx);
       current_user = final_decoded_users_again{current};
       root_duplicates = setdiff(find(decoded_roots==decoded_roots(current)), current);
       for compare_roots = root_duplicates
           if (~ismember(compare_roots, duplicate_set))
               compare_user = final_decoded_users_again{compare_roots};
               if (prod(compare_user==current_user))
                   duplicate_set = [duplicate_set, compare_roots];
               end
           end
       end       
   end
   
   if (~isempty(duplicate_set))
       final_decoded_users_again(duplicate_set) = [];
       num_decoded_users_again = size(final_decoded_users_again, 2);
   end

   
   num_decoded_users_total = num_decoded_users + num_decoded_users_again;

   for c_u_again = 1 : num_decoded_users_again
       decoded_msg = final_decoded_users_again{c_u_again};
       decoded_user = find(user_roots == decoded_msg(1));
       
       if (decoded_user)
           if (sum(prod(CRC_encoded_pos(:, decoded_user)==decoded_msg'))) % paths with the same roots are included
               num_correct_decoded = num_correct_decoded + 1;
           end
       end
       
   end

   num_missed_detection_again(trial) = num_active - num_correct_decoded;
   num_false_alarm_again(trial) = num_decoded_users_total - num_correct_decoded;

    
end

PUPE_MD = sum(num_missed_detection) / (num_trials*num_active);
PUPE_FA = sum(num_false_alarm) / (num_trials*num_active);

PUPE_MD_again = sum(num_missed_detection_again) / (num_trials*num_active);
PUPE_FA_again = sum(num_false_alarm_again) / (num_trials*num_active);

end

3.8 MAP_AMP_Hybrid_SPARCs_UMAC.m

function decoded_candidates = MAP_AMP_Hybrid_SPARCs_UMAC(y, n, M, K, extra_candidates,ordering, P, T)
%  this function corresponds to our inner decoder which consists of a
%  successive cancellation algorithm and a simplified AMP decoder.

% y : the received signal
% n : the number of channel uses, i.e., the length of transmitted codeword
% M: the section size of SPARCs
% K: the number of active users
% extra_candidates: the number of extra candidates
% ordering: chosen columns of the Hadamard matrix
% P: average power
% T: the maximum iterations for AMP decoding

allocated_P = sqrt(n*P);

decoded_candidates = zeros(1, K+extra_candidates);

n_2power = M;

% to avoid collisions among different active users, we use
% successive-cancellation-based MAP decoder as a preprocessing procedure.
K_sc = min(round(K/5), 10);
for l = 1 : K_sc
    beta_ = zeros(M, 1);
    inner_prod = FWHT_user_subsampling(2, y, n_2power, M, ordering);
    [~, decoded_info_pos_] = max(inner_prod);
    decoded_candidates(l) = decoded_info_pos_;
    beta_(decoded_info_pos_) = 1;
    y = y - sqrt(P)*FWHT_user_subsampling(1, beta_, n_2power, M, ordering);
end

     
% a simplified AMP decoder

K_remain = K - K_sc;
q = K_remain/M; %sparisity


z = y;
last_tao = 0;
beta = zeros(M, 1);

for t = 1 : T
    tao = sqrt(sum(z.^2)/n);
    if (abs(tao-last_tao) < 1e-6)
        break;
    end
    last_tao = tao;
                                        
    r = allocated_P*beta + FWHT_user_subsampling(2, z, n_2power, M, ordering)/sqrt(n);
     
    denoiser_1_part = exp(-(r-allocated_P).^2 / (2*tao^2)); 
    denoiser_0_part = exp(-r.^2/(2*tao^2));
    
    beta = q*denoiser_1_part ./ ((1-q)*denoiser_0_part+q*denoiser_1_part);
    z = y - sqrt(P)*FWHT_user_subsampling(1, beta, n_2power, M, ordering) ... 
                                       + (z/(tao.^2))*P*sum(beta.*(1-beta));
end

[~, topk_pos] = maxk(beta, K_remain+extra_candidates);
decoded_candidates(K_sc+1:K+extra_candidates) = topk_pos;

end

3.9 FWHT_user_subsampling.m

function x = FWHT_user_subsampling(mode, beta, n_2power, M, ordering)
% compute a general x = A*beta via a fast walsh-hadamard transform.

% ordering : the rows randomly chosen from Hadamard matrix.

switch (mode)
    
    % mode = 1 , for A*beta 
    case 1 
        beta_2power = zeros(n_2power, 1);
        beta_2power(n_2power - M + 1 : n_2power) = beta;
        xx = fwht_user(beta_2power')';
        x = xx(ordering);
        
    % mode = 2, for A'*z
    case 2
        beta_2power = zeros(n_2power, 1);
        beta_2power(ordering) = beta;
        xx = fwht_user(beta_2power')';
        x = xx(n_2power - M + 1 : n_2power);
end

end

3.10 fwht_user.c

/*=================================================================
 * Implement the fast walsh-hadamard transform using C to accelerate the speed. 
 * fwht_user.c
 * Copyright 2020-2020 Haiwen Cao.
 *	
 *=================================================================*/
#include <stdio.h>
#include <string.h>
#include "mex.h"


/* The computational routine */
void fwht_user(double *u, mwSize N)
{
    mwSize i, j, k, ij;
    double temp;

    for(i=N>>1; i>0; i>>=1) {
        for(k=0; k<N; k += 2*i) {
            for(j=k, ij=i+k; j<k+i; j++, ij++) {
                temp   = u[j];
                u[j]  += u[ij];
                u[ij] = temp - u[ij];
            }
        }
    }

}

3.11 fwht_user.mexw64

四、parameter configurations

4.1 configurations_UMAC_J_15_L_11_B_75 (Fig_3).mat

4.2 configurations_UMAC_J_14_L_7_B_50 (Fig_4).mat

4.3 configurations_UMAC_J_15_L_11_B_75 (Fig_5).mat

标签：Aided,decoded,Codes,end,Unsourced,current,num,round,roots
From： https://www.cnblogs.com/zhuxiaohuo/p/17923624.html

CRC-Aided Sparse Regression Codes for Unsourced Random Access