一、摘要

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

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_paths0);

remaining_roots_elimination(final_num_survival_paths0) = []; % 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 &gt; 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_rootsdecoded_roots(current)), current);

for compare_roots = root_duplicates

if (~ismember(compare_roots, duplicate_set))

compare_user = final_decoded_users{compare_roots};

if (prod(compare_usercurrent_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_trialsnum_active);

PUPE_FA = sum(num_false_alarm) / (num_trialsnum_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_rootsdecoded_roots(current)), current);

for compare_roots = root_duplicates

if (~ismember(compare_roots, duplicate_set))

compare_user = final_decoded_users{compare_roots};

if (prod(compare_usercurrent_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_rootsdecoded_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_usercurrent_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_trialsnum_active);

PUPE_FA = sum(num_false_alarm) / (num_trialsnum_active);
PUPE_MD_again = sum(num_missed_detection_again) / (num_trialsnum_active);

PUPE_FA_again = sum(num_false_alarm_again) / (num_trialsnum_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&gt;&gt;1; i&gt;0; i&gt;&gt;=1) {
    for(k=0; k&lt;N; k += 2*i) {
        for(j=k, ij=i+k; j&lt;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

本栏目推荐文章
• ACCESS 快速构建修改数据的窗体
• 解决前后端的跨域问题：Access to XMLHttpRequest at '**' from origin '**' has been blocked by CORS policy: No 'Access-Control-Allow-Origin' header is present on the requested resource.
• android编译kanzi 问题 (1) Caused by: java.io.FileNotFoundException: Error: Could not find or access Kanzi's Android libraries directory: C:\XXXXXXX
• No 'Access-Control-Allow-Origin' header is present on the requested resource'， 跨域访问的解决方法
• PV视角之3D检测模型Sparse4D系列
• CMU DLSys 课程笔记 2 - ML Refresher / Softmax Regression
• Access denied for user 'root'@'172.19.0.5' (using password: YES)...
• curl_easy_perform() failed: Problem with the SSL CA cert (path? access rights?)
• [ABC271G] Access Counter 题解
• AWS - Grant AWS EKS cluster access to Postgres and Redis using security group

Regression CRC-Aided Unsourced Access Sparseregression crc-aided unsourced access unsourced multiple activity access crc-aided unsourced sparse sparse4d sparse4 torch_sparse regression torch torch_scatter torch_sparse scatter