基于码分多址的多用户水下无线光通信（3）——解相关多用户检测

  继续上一篇博文，本文将介绍基于解相关的多用户检测算法。解相关检测器的优点是因不需要估计各个用户的接收信号幅值而具有抗远近效应的能力。常规的解相关检测器有运算量大和实时性差的缺点，本文针对异步CDMA的MAI主要来自干扰用户的相邻三个比特周期的特点，给出了基于相邻三个匹配滤波器输出数据的截断解相关检测算法。（我不知道怎么改公式里的字体，有的字母在公式中重复使用了，请根据上下文判断字母含义）

1、常规检测器

  假设已知各个用户的延时，且各个用户的延时满足$0\le {\tau }_1\le \cdots \le {\tau }_K<T_{\text{b}}$。无论是常规检测器还是多用户检测器，接收信号都要经过相关器进行解扩。在相关器中，待检测用户的扩频波形被重新生成并与接收信号进行相关运算。相关器可以通过匹配滤波技术实现，因此又被称为匹配滤波器。对匹配滤波器的输出的采样时刻与对应的待检测用户的信号延时同步，采样间隔为一个比特周期，该采样值是判决最有可能传输的信息比特的充分统计量。第$k$个用户的第$i$个比特的匹配滤波器输出采样值表示为 $$\begin{array}{rl}{y}_k[i]& =\sum _{n=0}^{L{N}_{\text{s}}-1}{s}_{\text{rx}}[n_0+(i-1)L{N}_{\text{s}}+q_k+n]s_k[n]\\ & =R_{k,k}(0)A_k{b}_k[i]+\sum _{m=-1}^1\sum _{\begin{array}{c}j=1\\ j\ne k\end{array}}^K{R}_{k,j}(m)A_j{b}_j[i-m]+v_k[i],\end{array}$$ 其中，$R_{k,k}(0)$是第$k$个用户的扩频波形在相对延时为$0$时的自相关值，$k\ne j$时的$R_{k,j}(m)$是两个不同用户的扩频波形之间的互相关值，$$R_{k,j}(m)=\sum _{n=-\infty }^{\infty }{s}_k[n-q_k]s_j[n+mL{N}_{\text{s}}-q_j],$$ $v_k[i]$表示匹配滤波器输出的噪声。$y_k[i]$公式中的第二行第一项表示有用的恢复数据，第二项表示匹配滤波器与其他用户的信号做相关运算产生的MAI。
  将针对各个用户的匹配滤波器的输出采样值按用户延时由短到长的顺序写入一个向量$\mathbf{y}(i)={\left[y_1[i],\cdots ,y_K[i]\right]}^{\text{T}}\in {\mathbb{R}}^{K\times 1}$，向量$\mathbf{y}(i)$表示为 $$\mathbf{y}(i)=\mathbf{R}(1)\mathbf{Q}\mathbf{b}(i-1)+\mathbf{R}(0)\mathbf{Q}\mathbf{b}(i)+\mathbf{R}(-1)\mathbf{Q}\mathbf{b}(i+1)+\mathbf{v}(i),$$ 其中，矩阵$\mathbf{R}(m)\in {\mathbb{R}}^{K\times K}$的第$(k,j)$个元素为$R_{k,j}(m)$，$\mathbf{R}(1)$是对角线为$0$的上三角矩阵，$\mathbf{R}(-1)$是对角线为$0$的下三角矩阵， b ( i ) = [ b 1 [ i ] , ⋯   , b K [ i ] ] T ∈ { + 1 , − 1 } K × 1 \boldsymbol{b}(i)=\left[b_1[i], \cdots,b_{K}[i]\right]^\text{T} \in \{+1,-1\}^{K\times 1} b(i)=[b1​[i],⋯,bK​[i]]T∈{+1,−1}K×1，$\mathbf{Q}=\text{diag}(A_1,\cdots ,A_K)\in {\mathbb{R}}^{K\times K}$，$\mathbf{v}(i)={\left[v_1[i],\cdots ,v_K[i]\right]}^{\text{T}}\in {\mathbb{R}}^{K\times 1}$。

  异步CDMA的常规检测器如上图所示，它由一组匹配滤波器和硬判决器构成，每个匹配滤波器对应一个用户。常规检测器直接对匹配滤波器的采样值做硬判决获得各个用户的信息比特，即 $$\hat{\mathbf{b}}(i)=\text{sgn}\left(\mathbf{y}(i)\right).$$ 常规检测器遵循单用户检测策略，它的每个分支只检测一个用户，并且把其他干扰用户的信号视为噪声。因此，常规检测器没有利用多用户的信息或进行联合信号处理。当系统中用户数量多且扩频码非正交时，会有严重的MAI，这对常规检测器的检测性能产生严重负面影响。此外，当各个用户的信号以不同功率到达接收机，即存在远近效应时，信号强的用户会加剧信号弱的用户的MAI，弱的信号会被强的信号淹没。

function [user_bits,value] = conventional_mult_demod(rec_data,ss_code,sf,L_bit,K,b,sps,delay,isShape,isDsShape)
%Conventional single-user detector
% rec_data 接收信号
% ss_code 扩频码
% sf 扩频因子
% L_bit 信息比特数
% b 成型滤波抽头系数
% K 用户数
% sps 上采样率
% delay 传输延时
% isShape 发射信号是否成型滤波 1 滤波, 0 不滤波
% isDsShape 用于解扩的本地码是否成型滤波 1 滤波, 0 不滤波
if nargin < 10
    isDsShape = 0;
end
if isDsShape == 1 && isShape == 1
    ss_code_rcos = upfirdn(ss_code',b,sps)';
    ss_code_rcos = ss_code_rcos(1:K,round(length(b)/2)-round(sps/2)+1:round(length(b)/2)-round(sps/2)+sps*sf);
end
value = zeros(K,L_bit);
for k = 1:K
    temp = reshape(rec_data(delay(k):L_bit*sps*sf+delay(k)-1)',[],L_bit)';
    if isDsShape == 1 && isShape == 1
        value(k,:) = (temp*ss_code_rcos(k,:)')';
    else
        value(k,:) = (temp*rectpulse(ss_code(k,:)',sps))';
    end 
end
user_bits = sign(value);
user_bits(user_bits == 0) = 1;
user_bits = (user_bits+1)/2;
end

2、基于解相关的多用户检测器

  解相关检测器是在常规检测器的匹配滤波器组后面加了一级线性变换，以消除各个用户扩频波形之间的相关性。
  常规的解相关检测器是在获取多个比特周期的数据后再检测每个用户的信息比特，假设解相关检测器每$N$个比特周期做一次解相关检测。将匹配滤波器组的$N$个采样值写入一个列向量${\mathbf{y}}_N(i)\in {\mathbb{R}}^{NK\times 1}$，${\mathbf{y}}_N(i)$表示为 $${\mathbf{y}}_N(i)=\left[\begin{array}{c}\mathbf{y}(i)\\ \mathbf{y}(i+1)\\ ⋮\\ \mathbf{y}(i+N-1)\end{array}\right].$$ ${\mathbf{y}}_N$的矩阵表达为 $${\mathbf{y}}_N(i)={\mathbf{R}}_N{\mathbf{A}}_N{\mathbf{b}}_N(i)+{\mathbf{v}}_N(i),$$ 其中，${\mathbf{R}}_N\in {\mathbb{R}}^{NK\times NK}$是一个分块Toeplitz矩阵 $${\mathbf{R}}_N=\left[\begin{array}{ccccc}\mathbf{R}(0)& \mathbf{R}(-1)& {\text{0}}_{K\times K}& \cdots & {\text{0}}_{K\times K}\\ \mathbf{R}(1)& \mathbf{R}(0)& \mathbf{R}(-1)& \ddots & ⋮\\ {\text{0}}_{K\times K}& \mathbf{R}(1)& \mathbf{R}(0)& \ddots & {\text{0}}_{K\times K}\\ ⋮& \ddots & \ddots & \ddots & \mathbf{R}(-1)\\ {\text{0}}_{K\times K}& \cdots & {\text{0}}_{K\times K}& \mathbf{R}(1)& \mathbf{R}(0)\end{array}\right],$$ ${\mathbf{A}}_N\in {\mathbb{R}}^{NK\times NK}$表示由信号幅值构成的对角矩阵 $${\mathbf{A}}_N=\left[\begin{array}{ccc}\mathbf{Q}& & \\ & \ddots & \\ & & \mathbf{Q}\end{array}\right],$$ b N ( i ) ∈ { + 1 , − 1 } N K × 1 \boldsymbol{b}_N(i)\in \{+1,-1\}^{NK\times 1} bN​(i)∈{+1,−1}NK×1 包含$K$个用户的共$NK$个信息比特 $${\mathbf{b}}_N(i)=\left[\begin{array}{c}\mathbf{b}(i)\\ \mathbf{b}(i+1)\\ ⋮\\ \mathbf{b}(i+N-1)\end{array}\right],$$ ${\mathbf{v}}_N(i)$包含匹配滤波器输出的噪声以及来自$\mathbf{b}(i-1)$和$\mathbf{b}(i+N)$的干扰。解相关检测器的输出结果表示为 $$\begin{array}{rl}{\hat{\mathbf{b}}}_N(i)& =\text{sgn}\left({\mathbf{R}}_N^{-1}{\mathbf{y}}_N(i)\right)\\ & =\text{sgn}\left({\mathbf{A}}_N{\mathbf{b}}_N(i)+{\mathbf{R}}_N^{-1}{\mathbf{v}}_N(i)\right).\end{array}$$ 由上式可以看出，解相关检测器不需要估计用户的信号的幅值，因此解相关检测器适用于信号动态范围较大的水下环境。
  如果$N$非常大，则求${\mathbf{R}}_N^{-1}$需要大量的计算，这将使得解相关检测器的运算量和实时性难以被接受。根据前面的介绍可知，异步CDMA的MAI主要来自干扰用户的相邻三个比特周期。可以截取出每个用户的三个相邻比特对应的匹配滤波器输出作为一组用于解相关检测，同时为了保证检测感兴趣的信息比特所需的信息包含在接收到的信号中，判决结果只取中间比特，下图给出了这种截断解相关检测器的示意图。

  将每个用户的三个相邻比特对应的匹配滤波器输出写入一个向量$\mathbf{\xi }(i)\in {\mathbb{R}}^{3K\times 1}$，$\mathbf{\xi }(i)$表示为 $$\begin{array}{rl}\mathbf{\xi }(i)& =\left[\begin{array}{c}\mathbf{y}(i-1)\\ \mathbf{y}(i)\\ \mathbf{y}(i+1)\end{array}\right]\\ & =\mathbf{R}\mathbf{A}\mathbf{x}(i)+\mathbf{\nu }(i),\end{array}$$其中，$\mathbf{R}\in {\mathbb{R}}^{3K\times 3K}$为 $$\mathbf{R}=\left[\begin{array}{ccc}\mathbf{R}(0)& \mathbf{R}(-1)& {\text{0}}_{K\times K}\\ \mathbf{R}(1)& \mathbf{R}(0)& \mathbf{R}(-1)\\ {\text{0}}_{K\times K}& \mathbf{R}(1)& \mathbf{R}(0)\end{array}\right],$$ $\mathbf{A}\in {\mathbb{R}}^{3K\times 3K}$表示为 $$\mathbf{A}=\left[\begin{array}{ccc}\mathbf{Q}& & \\ & \mathbf{Q}& \\ & & \mathbf{Q}\end{array}\right],$$ x ( i ) ∈ { + 1 , − 1 } 3 K × 1 \boldsymbol{x}(i)\in \{+1,-1\}^{3K\times 1} x(i)∈{+1,−1}3K×1 包含$K$个用户的共$NK$个信息比特$$\mathbf{x}(i)=\left[\begin{array}{c}\mathbf{b}(i-1)\\ \mathbf{b}(i)\\ \mathbf{b}(i+1)\end{array}\right],$$ $\mathbf{\nu }(i)$包含匹配滤波器输出的噪声以及来自$\mathbf{b}(i-2)$和$\mathbf{b}(i+2)$的干扰 $$\mathbf{\nu }(i)=\left[\begin{array}{c}\mathbf{R}(1)\mathbf{Q}\mathbf{b}(i-2)+\mathbf{v}(i-1)\\ \mathbf{v}(i)\\ \mathbf{R}(-1)\mathbf{Q}\mathbf{b}(i+2)+\mathbf{v}(i+1)\end{array}\right].$$ 对于每个用户的第$i$个信息比特，截断解相关检测器的输出向量表示为 $$\begin{array}{rl}\hat{\mathbf{b}}(i)& =\text{sgn}\left(\mathbf{S}{\mathbf{R}}^{-1}\mathbf{\xi }(i)\right)\\ & =\text{sgn}\left(\mathbf{S}\left(\mathbf{A}\mathbf{x}(i)+{\mathbf{R}}^{-1}\mathbf{\nu }(i)\right)\right),\end{array}$$ 其中，$\mathbf{S}=\left[{\text{0}}_{K\times K},{\mathbf{I}}_{K\times K},{\text{0}}_{K\times K}\right]$是一个选择矩阵，用于选出每个用户的第$i$个信息比特，${\mathbf{I}}_{K\times K}$是$K$维单位矩阵。

function [user_bits,value] = decorrelating_mult_demod(rec_data,ss_code,sf,L_bit,K,b,sps,delay,isShape,isDsShape)
%Decorrelating multiuser detector
% rec_data 接收信号
% ss_code 扩频码
% sf 扩频因子
% L_bit 信息比特数
% K 用户数
% b 成型滤波抽头系数
% sps 上采样率
% delay 传输延时
% isShape 发射信号是否成型滤波 1 滤波, 0 不滤波
% isDsShape 用于解扩的本地码是否成型滤波 1 滤波, 0 不滤波
if nargin < 10
    isDsShape = 0;
end
ss_code_rect = rectpulse(ss_code',sps)';
% ss_code_rect = ss_code_rect./vecnorm(ss_code_rect,2,2); 
if isShape == 1
    ss_code_rcos = upfirdn(ss_code',b,sps)';
    ss_code_rcos = ss_code_rcos(:,round(length(b)/2)-round(sps/2)+1:round(length(b)/2)-round(sps/2)+sps*sf);
end
% 计算扩频码的互相关矩阵 R_{11}, R_{12}, ..., R_{1K}, R_{21}, ..., R_{KK},第一个下标对应根升余弦滤波的波形，第二个下标对应矩形成型波形
R_ss_code = zeros(2*length(ss_code_rect)-1,K*K); 
for k = 1:K
    for l = 1:K
        if isDsShape == 1 && isShape == 1
            R_ss_code(:,(k-1)*K+l) = xcorr(ss_code_rcos(k,:)',ss_code_rcos(l,:)','none');
        elseif isDsShape == 0 && isShape == 1
            R_ss_code(:,(k-1)*K+l) = xcorr(ss_code_rcos(k,:)',ss_code_rect(l,:)','none');
        else
            R_ss_code(:,(k-1)*K+l) = xcorr(ss_code_rect(k,:)',ss_code_rect(l,:)','none');
        end
    end
end
[~,I] = sort(delay,'ascend');  %到达时间排序
% 构造解相关矩阵
R_1 = zeros(K); %R(-1) lower triangular matrix 
R0 = zeros(K); %R(0)
R1 = zeros(K); %R(1) upper triangular matrix 
for k = 1:K
    k_ = I(k); %目标用户
    for l = 1:K
        l_ = I(l); % l_~=k_为干扰用户
        tao = delay(k_)-delay(l_);
        if tao == sps*sf || tao == -sps*sf
            tao = 0;
        end
        R0(k,l) = R_ss_code(sf*sps+tao,(l_-1)*K+k_);
        if l < k && (sf*sps+(tao-sf*sps) > 0 && sf*sps+(tao-sf*sps) < 2*sf*sps)
            R_1(k,l) = R_ss_code(sf*sps+(tao-sf*sps),(l_-1)*K+k_);
        end
        if l > k && (sf*sps+(sf*sps+tao) > 0 && sf*sps+(sf*sps+tao) < 2*sf*sps)
            R1(k,l) = R_ss_code(sf*sps+(sf*sps+tao),(l_-1)*K+k_);
        end
    end
end
R = [R0,      R_1, zeros(K);
     R1,      R0,  R_1;
     zeros(K),R1,  R0];
R_ = pinv(R);
% 计算每个用户各个bit对应的相关器输出，按到达时间早晚排列
y = zeros(K,L_bit);
for k = 1:K
    ind = I(k);
    temp = reshape(rec_data(delay(ind):L_bit*sps*sf+delay(ind)-1)',[],L_bit)';
    if isDsShape == 1 && isShape == 1
        y(k,:) = (temp*ss_code_rcos(ind,:)')';
    else
        y(k,:) = (temp*rectpulse(ss_code(ind,:)',sps))';
    end 
end
y = reshape(y,[],1);
value_temp = zeros(K*L_bit,1);
R_ = R_(K+1:2*K,:);
for n = 1:L_bit
    if n == 1
        temp = R_*[zeros(K,1);y((n-1)*K+1:(n+1)*K)];
    elseif n == L_bit
        temp = R_*[y((n-2)*K+1:n*K);zeros(K,1)];
    else
        temp = R_*y((n-2)*K+1:(n+1)*K);
    end
    value_temp((n-1)*K+1:n*K) = temp;
end
value = zeros(K,L_bit);
value_temp = reshape(value_temp,K,[]);
value(I,:) = value_temp;
user_bits = sign(value);
user_bits(user_bits == 0) = 1;
user_bits = (user_bits+1)/2;
end

3、算法仿真

  下面给一个仿真的顶层代码，遍历参数有信噪比和信干比，感兴趣的读者可以试一下看看效果。

date = '5_28_';
if(~exist(['.\',date,'sim_data'],'dir'))
    mkdir(['.\',date,'sim_data']);
end
K = 3; % 用户数
Ns = 10 ; % samples/chip
isShape = 1; %是否成型滤波 1 滤波, 0 不滤波
isDC = 0; % 接收机直流耦合 1 直流耦合, 0 交流耦合, 可以直流耦合接收，后面在代码里去直
isDsShape = 0; % 解扩时的本地码是否成型滤波 1 滤波, 0 不滤波
isEst = 1; % 是否信道估计
isTest = 0; % 测试
Target_User = 1; % 目标用户
noise_power = 22:-2:8; % noise power dBW
% noise_power = 10; % noise power dBW
target_user_power = 0; % AC power (variance) dBW
interference_user_power = [-10,0,10,20]; % AC power (variance) dBW
M = 200; % 快拍数
% 扩频码的PN序列多项式和初始值
ss_polynomial = [1 0 1 0 0 1;    % z^5+z^3+1
                 1 1 1 1 0 1;    % z^5+z^4+z^3+z^2+1
                 1 1 0 1 1 1];   % z^5+z^4+z^2+z^1+1
ss_init_state = [1 0 1 0 1;
                 1 0 1 0 1;
                 1 0 1 0 1];
if isShape == 1
    shape = '_rcos';
else
    shape = [];
end
if isDsShape == 1 && isShape == 1
    ds_shape = '_rcos';
else
    ds_shape = [];
end
if isEst == 1
    est = '_est';
else
    est = [];
end
if isTest == 1
    test = '_test1';
else
    test = [];
end
% 用户发送数据bit的PN序列多项式和初始值
bit_polynomial = [1,zeros(1,16),1,0,0,1; % z^20+z^3+1
                  1,zeros(1,10),1,zeros(1,3),1,0,1,0,0,1; % z^20+z^9+z^5+z^3+1
                  1,1,zeros(1,14),1,1,0,0,1];% z^20+z^19+z^4+z^3+1
bit_init_state = [repmat([1,0],1,10);
                  repmat([1,0],1,10);
                  repmat([1,0],1,10)]; 
% 用户数据bit帧头多项式和初始值
head_polynomial = [1 0 0 0 0 0 1 1;     % z^7+z+1
                   1 0 0 0 1 0 0 1;     % z^7+z^3+1
                   1 0 0 0 1 1 1 0];    % z^7+z^3+z^2+z+1
head_init_state = [1 0 1 0 1 0 1;
                   1 0 1 0 1 0 1;
                   1 0 1 0 1 0 1];              
L_ss = 2^(length(ss_polynomial)-1)-1; % length of spread spectrum pn sequence, spreading factor
L_head = 2^(length(head_polynomial)-1)-1;
L_bits = 1e5;
Times = 5;
delay_array = 0:8:L_ss*Ns-1;
%% 生成发送数据
ss_code = zeros(K,L_ss);
user_bits = zeros(K,L_bits);
user_head = zeros(K,L_head);
user_ss_data = zeros(K,(L_bits+L_head)*L_ss);
for k = 1:K
    % 扩频码
    ss_code(k,:) = 2*PnCode(ss_polynomial(k,:),ss_init_state(k,:))-1;
    % 用户数据bit
    user_bits_temp = 2*PnCode(bit_polynomial(k,:),bit_init_state(k,:))-1;
    user_bits(k,:) = user_bits_temp(1:L_bits);
    % 帧头
    user_head(k,:) = 2*PnCode(head_polynomial(k,:),head_init_state(k,:))-1;
    user_data_upsample = rectpulse([user_head(k,:),user_bits(k,:)],L_ss);
    user_ss_data(k,:) = user_data_upsample.*repmat(ss_code(k,:),1,L_bits+L_head); 
end
% 上采样，成型滤波
if isShape == 1
    sps = Ns; % upsample rate
    span = 6;
    rolloff = 0.5;
    b = rcosdesign(rolloff,span,sps,'sqrt');% 设计根升余弦滤波器
    % 成型滤波
    user_ss_data = upfirdn(user_ss_data',b,sps)';
else
    b = 1;
    sps = Ns;
    user_ss_data = rectpulse(user_ss_data',sps)';% 矩形成型
end
%归一化
user_ss_data = user_ss_data./vecnorm(user_ss_data,2,2).*sqrt(length(user_ss_data)); %功率归一化
clear user_bits_temp;
clear user_bits_upsample;
%% 接收
L_data = length(user_ss_data);
L_sample = L_data+Ns*L_ss;
BER = zeros(1+2*length(interference_user_power),length(noise_power)); % 记录一个单用户和多用户时的三种方法的误码率
user_ss_data(Target_User,:) = user_ss_data(Target_User,:)*sqrt(10^(target_user_power/10));
for t = 1:Times
    for n = 1:length(noise_power)  
        background_noise = wgn(1,L_sample,noise_power(n));% background noise
        single_user_rec_data = user_ss_data(Target_User,:)+background_noise(1:L_data);
        single_user_rec_data = single_user_rec_data-mean(single_user_rec_data); % DC block
        if isShape == 1
            delay = round(length(b)/2)-round(sps/2)+1; 
        else
            delay = 1; 
        end
        [single_user_conv_demod_bits,~] = conventional_mult_demod(single_user_rec_data,ss_code,L_ss,L_bits+L_head,1,b,Ns,delay,isShape,isDsShape);
        [~,ber] = biterr(single_user_conv_demod_bits(L_head+1:end),(user_bits(Target_User,:)+1)/2,[],'row-wise');
        BER(1,n) = BER(1,n)+ber;
        for d = 1:length(delay_array)
            send_data = zeros(K,L_sample);
            send_data(Target_User,delay_array(d)+1:delay_array(d)+L_data) = user_ss_data(Target_User,:);
            for p = 1:length(interference_user_power)
                background_noise = wgn(1,L_sample,noise_power(n));% background noise
                for k = 1:K
                    %模拟发送信号延时
                    if k~=Target_User
                        send_data(k,1:L_data) = user_ss_data(k,:).*sqrt(10^(interference_user_power(p)/10));
                    end 
                end
%                 near_far_ratio(p) = 10*log10(var(send_data(2,:))/var(send_data(Target_User,:)));
                rec_data = sum(send_data,1)+background_noise;
                rec_data = rec_data-mean(rec_data); % DC block
                if isEst == 1
                    % channel estimation
                    [delay,~] = subspace_geo_channel_estimation(rec_data,ss_code,L_ss,M,L_bits,K,b,Ns,isShape,1);
                    ind  = find(delay < length(b)/2-round(sps/2)-Ns/2);
                    if (isShape == 1) && (~isempty(ind))
                        delay(ind) = delay(ind)+Ns*L_ss;
                    end
                    delay = delay+1;
                else
                   %假设已经准确同步
                    if isShape == 1
                        delay = round(length(b)/2)-round(sps/2)+[delay_array(d),0,0]+1;
                    else
                        delay = [delay_array(d),0,0]+1;
                    end
                end
                % 常规检测器
                [conv_demod_bits,conv_demod_value] = conventional_mult_demod(rec_data,ss_code,L_ss,L_bits+L_head,K,b,Ns,delay,isShape,isDsShape);
                [~,conv_ber] = biterr(conv_demod_bits(:,L_head+1:end),(user_bits+1)/2,[],'row-wise');
                BER(1+(p-1)*2+1,n) = BER(1+(p-1)*2+1,n)+conv_ber(Target_User);
                % 解相关检测器
                [decorr_demod_bits,decorr_demod_value] = decorrelating_mult_demod(rec_data,ss_code,L_ss,L_bits+L_head,K,b,Ns,delay,isShape,isDsShape);
                [~,decorr_ber] = biterr(decorr_demod_bits(:,L_head+1:end),(user_bits+1)/2,[],'row-wise');
                BER(1+(p-1)*2+2,n) = BER(1+(p-1)*2+2,n)+decorr_ber(Target_User);
            end
        end
    end
end
SNR = target_user_power-noise_power;
EbN0 = SNR-10*log10(1/L_ss)+10*log10(Ns/2);
BER = BER./t;
BER(2:end,:) = BER(2:end,:)./d;
MAI = interference_user_power-target_user_power; % near far ratio
save(['.\',date,'sim_data\',date,'sim_ber_avr',est,test],'BER','SNR','EbN0');

  代码有点多，可能有的函数没贴上来，缺代码的话请留言、私信或者点此下载未删减的全套代码。
  代码中设置的三个m序列互为优选对，互为优选对的m序列具有三值互相关函数，且互相关函数值比较小，这对常规检测器来说是有利的。

  上图给出了不同MAI条件下，常规检测器（标记为“Conv.”）和解相关检测器（标记为“Decor.”）的BER与$E_{\text{b}}/N_0$关系关系的仿真结果。单用户（标记为“Single User”）时采用的是常规检测器，多用户检测器的BER曲线越接近单用户的BER曲线则说明多用户检测器抗MAI和远近效应的性能越好。常规检测器的性能对各个用户的相对功率有很强的依赖性，常规检测器的性能随着MAI的增强而迅速恶化，MAI超过$10$ dB时，常规检测器的性能实际上受到MAI限制。解相关检测器的误码性能受MAI影响小，接近单用户的BER。

参考文献

[1] LUPAS R, VERDU S. Linear multiuser detectors for synchronous code-division multipleaccess channels[J]. IEEE Transactions on Information Theory, 1989, 35(1): 123-136.

[2] LUPAS R, VERDU S. Near-far resistance of multiuser detectors in asynchronous channels[J]. IEEE Transactions on Communications, 1990, 38(4): 496-508.