基于码分多址的多用户水下无线光通信（2）——系统模型和基于子空间的延时估计

  本文首先介绍了基于CDMA的多用户UOWC系统模型，并给出了多用户收发信号的数学模型。然后介绍基于子空间的延时估计算法，该算法只需要已知所有用户的扩频码，然后根据扩频波形的循环移位在观测空间的信号子空间上的投影进行延时估计。

1、基于CDMA的多用户UOWC系统模型

  首先介绍基于CDMA的多用户UOWC系统模型，系统框图如下图所示。

  该系统包括发送端、UOWC信道和接收端。该系统包含$K$个用户，每个用户配备一个LED光源，且为其分配了独一无二的m序列作为扩频码。与射频通信不同的是，LED属于非相干光源，因此发送端只能采用强度调制（调制LED发光的亮灭或者亮的程度），并且LED只能发送单极性的实信号，需要用Bias-T将双极性信号叠加上直流偏置以使LED工作在输出光功率和输入电压关系曲线的线性区间。在接收端采用隔直的雪崩光电二极管（Avalanche Photodiode，APD）作为接入点（Access Point，AP）的光电探测器。光信号被隔直APD转换为电信号并去除直流分量，交流电信号经ADC采样后获得离散数字信号。
  假设这$K$个用户有相同的比特周期$T_{\text{b}}$和码片周期$T_{\text{c}}$。比特周期与码片周期的关系为$$T_{\text{c}}=\frac{T_{\text{b}}}{L},$$ 其中$L$为扩频因子。第$k$个用户发送的信号可以表示为$$s_{\text{tx},k}(t)=\sum _{i=-\infty }^{\infty }{P}_k{b}_k[i]s_k(t-i{T}_{\text{b}})+m_k,$$其中，$i$是信息比特的索引，$P_k$表示第$k$个用户发送信号的幅度， b k [ i ] ∈ { + 1 , − 1 } b_k [i]\in\{+1,-1\} bk​[i]∈{+1,−1}表示第$k$个用户发送的第$i$个比特，$s_k(t)$表示第$k$个用户扩频波形，它在区间$[0,T_{\text{b}}]$以外的值为$0$，$m_k$是加在第$k$个用户的LED上的直流偏置。
  所有用户的光信号经过水下信道后在接收端叠加，叠加的光信号经过隔直APD的光电转换和去直后变成双极性电信号。ADC以$T_{\text{s}}=T_{\text{c}}/N_{\text{s}}$为采样间隔对电信号采样获得离散数字信号，其中$N_{\text{s}}$是上采样率。$n{T}_{\text{s}}$时刻的采样信号表示为$$s_{\text{rx}}[n]=\sum _{k=1}^K\sum _{i=-\infty }^{\infty }{A}_k{b}_k[i]s_k\left[n-iL{N}_{\text{s}}-q_k\right]+w[n],$$ 其中，$A_k=P_k{h}_k$表示电信号的幅值，$h_k$表示包括电光转换、水下信道功率损耗和光电转换的等效信道增益（不考虑多径传输），$s_k[n]=s_k(n{T}_{\text{s}})$表示扩频波形的采样值，$q_k$表示延时对应的采样点数量，$w[n]$是均值为$0$方差为${\sigma }^2$的高斯白噪声。

2、延时估计

  基于异步CDMA的多用户UOWC系统的信号如下图所示。

  假设数据以连续比特进行传输，每个用户的扩频码已知。由于各个用户和接收机之间是异步的，每个用户的信号在不同时刻到达接收机，接收机一开始不知道每个用户的比特的确切位置。延时估计的目标是确定每个用户的信号解扩时积分区间的起始点，也可以视为位同步。以$n_0{T}_{\text{s}}$为参考采样时刻，对于第$k$个用户，第$n_0$个采样点所处的比特记为$b_k[0]$，待估计的延时就是从第$n_0$个采样点到$b_k[1]$比特第一个采样点之间的时间间隔。从上图中容易看出，待估计的延时不超过一个比特周期。即使某个用户的延时${\tau }_k$超过了一个比特周期，也可以对${\tau }_k$按比特周期取模将其限制在一个比特周期以内，即$({\tau }_k\text{mod}{T}_{\text{b}})\in [0,T_{\text{b}})$。因此，可以不失一般性的假设每个用户的延时不超过一个比特周期，即${\tau }_k\in [0,T_{\text{b}})$。由于采样间隔$T_{\text{s}}$是系统中最小的时间单元，估计延时${\tau }_k$就是估计${\tau }_k$对应的采样点数量$q_k$：$$q_k=⟨\frac{{\tau }_k}{T_{\text{s}}}⟩\text{mod}(L{N}_{\text{s}}),$$其中，$⟨\cdot ⟩$表示四舍五入取整，$q_k$从集合 { 0 , 1 , ⋯   , L N s − 1 } \{0,1,\cdots,LN_\text{s}-1\} {0,1,⋯,LNs​−1}中取值。

2.1 滑动相关法延时估计

  常规的延时估计方法是滑动相关法，该方法用滑动相关器与接收信号做相关运算，通过最大化滑动窗口内的相关值来估计信号延时。滑动相关器是一个与待估计用户的扩频波形匹配的FIR结构滤波器，第$k$个用户的滑动相关器的滤波器系数就是其扩频波形在区间$[0,T_{\text{b}}]$内的采样值，用向量表示为 $${\mathbf{s}}_k=\left[s_k[0],s_k[1],\cdots ,s_k[L{N}_{\text{s}}-1]\right].$$ 假设用持续$M$个比特周期的接收信号估计第$k$个用户的延时，滑动相关法的公式表示为 q ^ k = arg ⁡ max ⁡ q ∈ { 0 , ⋯   , L N s − 1 } ∑ i = 0 M − 1 ∣ s k r ( q , i ) ∣ , \hat{q}_k = \mathop{\arg\max}\limits_{q \in \{0,\cdots,LN_\text{s}-1\}}\sum_{i=0}^{M-1}{\left|\boldsymbol{s}_k\boldsymbol{r}(q,i)\right|}, q^​k​=q∈{0,⋯,LNs​−1}argmax​i=0∑M−1​∣sk​r(q,i)∣, 其中，$$\mathbf{r}(q,i)=\left[\begin{array}{c}{s}_{\text{rx}}[n_0+q+iL{N}_{\text{s}}]\\ s_{\text{rx}}[n_0+q+iL{N}_{\text{s}}+1]\\ ⋮\\ s_{\text{rx}}[n_0+q+(i+1)L{N}_{\text{s}}-1]\end{array}\right]\in {\mathbb{R}}^{L{N}_{\text{s}}\times 1}.$$
  滑动相关法将MAI视为噪声，它能够在单用户或者多用户扩频波形相互正交的情况下工作。然而，在有MAI，尤其是存在远近效应的情况下，互相关峰可能大于自相关峰，导致滑动相关法不能正常工作。不准确的延时估计结果将使多用户检测器的检测性能变差，因此有必要研究抗远近效应的延时估计算法。

function [delay,max_corr] = corr_channel_estimation(rec_data,ss_code,sf,M,L_bit,K,b,sps,isShape)
% rec_data 接收信号
% ss_code 扩频码
% sf 扩频因子
% M 快拍数
% L_bit 信息比特数
% K 用户数
% b 成型滤波抽头系数
% sps 上采样率
% isShape 发射信号是否成型滤波 1 滤波, 0 不滤波
if M >= L_bit
    M = L_bit-1;
end
if isShape == 1
    ss_code_upsample = upfirdn(ss_code',b,sps)';
    ss_code_upsample = ss_code_upsample(:,round(length(b)/2)-round(sps/2)+1:round(length(b)/2)-round(sps/2)+sps*sf);
else
    ss_code_upsample = rectpulse(ss_code,sps);
end
delay = zeros(K,1);
max_corr = zeros(K,1);
corr = zeros(1,M*sps*sf);
for k = 1:K
    for n = 1:M*sps*sf
        corr(n) = rec_data(n:n+sps*sf-1)*ss_code_upsample(k,:)';
    end
    corr_temp = reshape(abs(corr'),sps*sf,[]);
    [max_corr(k),delay(k)] = max(mean(corr_temp,2));
end
delay = delay-1;
end

2.2 基于子空间的延时估计

  这个方法和阵列信号处理中MUSIC算法很像，阵列信号处理中是用多个天线接获得正弦信号的不同相位，而这里是一个光电探测器长时间采样获得扩频序列的不同相位。
  令一个观测窗口的宽度等于一个比特周期$T_{\text{b}}$，有$L{N}_{\text{s}}$个采样点。从第$n_0$个采样点开始，将每$L{N}_{\text{s}}$个采样点写入一个观测向量$\mathbf{z}(i)\in {\mathbb{R}}^{L{N}_{\text{s}}\times 1}$，$\mathbf{z}(i)$表示为 $$\mathbf{z}(i)=\left[\begin{array}{c}{s}_{\text{rx}}[n_0+iL{N}_{\text{s}}]\\ s_{\text{rx}}[n_0+iL{N}_{\text{s}}+1]\\ ⋮\\ s_{\text{rx}}[n_0+(i+1)L{N}_{\text{s}}-1]\end{array}\right],$$ 其中，$i$表示观测窗口起始位置在第$i$个比特。虽然每个观测窗口对应一个比特周期，但是$\mathbf{z}(i)$是在不考虑比特同步的情况下获得的，每个观测窗口可能包含每个活跃用户的当前比特的末尾和下一个比特的开头。根据${\mathbf{s}}_k$和延时采样点数量$q_k$可以定义两个向量${\mathbf{u}}_k^{\text{r}}\in {\mathbb{R}}^{L{N}_{\text{s}}\times 1}$和${\mathbf{u}}_k^{\text{l}}\in {\mathbb{R}}^{L{N}_{\text{s}}\times 1}$，其中${\mathbf{u}}_k^{\text{r}}$由${\mathbf{s}}_k$的右半部分和补0组成，${\mathbf{u}}_k^{\text{r}}$表示为 $${\mathbf{u}}_k^{\text{r}}=\{\begin{array}{ll}{\text{0}}_{L{N}_{\text{s}}\times 1},& q_k=0\\ {\left[s[L{N}_{\text{s}}-q_k],\cdots ,s[L{N}_{\text{s}}-1],0,\cdots ,0\right]}^{\text{T}},& 1\le q_k\le L{N}_{\text{s}}-1\end{array},$$ ${\mathbf{u}}_k^{\text{l}}$由补0和${\mathbf{s}}_k$的左半部分组成，${\mathbf{u}}_k^{\text{l}}$表示为 $${\mathbf{u}}_k^{\text{l}}=\{\begin{array}{ll}{\mathbf{s}}_k^{\text{T}},& q_k=0\\ {\left[0,\cdots ,0,s[0],\cdots ,s[L{N}_{\text{s}}-1-q_k]\right]}^{\text{T}},& 1\le q_k\le L{N}_{\text{s}}-1\end{array}.$$ 将观测向量$\mathbf{z}(i)$进一步表示为 $$\begin{array}{rl}\mathbf{z}(i)& =\sum _{k=1}^K\left[{\mathbf{u}}_k^{\text{r}}{A}_k{b}_k[i]+{\mathbf{u}}_k^{\text{l}}{A}_k{b}_k[i+1]\right]+\mathbf{w}(i)\\ & =\mathbf{U}\mathbf{A}\mathbf{d}(i)+\mathbf{w}(i),\end{array}$$ 其中，$\mathbf{U}=[{\mathbf{u}}_1^{\text{r}},{\mathbf{u}}_1^{\text{l}},\cdots ,{\mathbf{u}}_K^{\text{r}},{\mathbf{u}}_K^{\text{l}}]\in {\mathbb{R}}^{L{N}_{\text{s}}\times 2K}$，$\mathbf{A}=\text{diag}(A_1,A_1,\cdots ,A_K,A_K)\in {\mathbb{R}}^{2K\times 2K}$， d ( i ) = [ b 1 [ i ] , b 1 [ i + 1 ] , ⋯   , b K [ i ] , b K [ i + 1 ] ] T ∈ { + 1 , − 1 } 2 K × 1 \boldsymbol{d}(i) = \left[b_1[i] , b_1[i+1], \cdots , b_{K}[i] , b_{K}[i+1] \right]^\text{T} \in \{+1,-1\}^{2K\times 1} d(i)=[b1​[i],b1​[i+1],⋯,bK​[i],bK​[i+1]]T∈{+1,−1}2K×1，$\mathbf{w}(i)={\left[w[n_0+iL{N}_{\text{s}}],\cdots ,w[n_0+(i+1)L{N}_{\text{s}}-1]\right]}^{\text{T}}\in {\mathbb{R}}^{L{N}_{\text{s}}\times 1}$是高斯白噪声向量。延时估计是利用观测向量$\mathbf{z}(i)$识别出发送信号的用户的扩频码并估计信号延时$q_k$。
  假设每个用户发送独立同分布的信息比特，噪声与发送的信息无关，并且待估计参数在估计阶段保持不变。首先用多组观测向量$\mathbf{z}(i)$估计自相关矩阵 $${\hat{\mathbf{R}}}_{\text{z}}=\frac{1}{M}\sum _{i=0}^{M-1}\mathbf{z}(i){\mathbf{z}}^{\text{T}}(i),$$ 其中，$M$是观测窗口数量，即快拍数。下一步对${\hat{\mathbf{R}}}_{\text{z}}$做特征值分解，并根据特征值大小将特征向量分成信号子空间和噪声子空间，表示为 $$\begin{array}{rl}{\hat{\mathbf{R}}}_{\text{z}}& =[{\hat{\mathbf{V}}}_{\text{s}},{\hat{\mathbf{V}}}_{\text{n}}]\left[\begin{array}{cc}{\hat{\mathbf{\Lambda }}}_{\text{s}}& \text{0}\\ \text{0}& {\hat{\mathbf{\Lambda }}}_{\text{n}}\end{array}\right]\left[\begin{array}{c}{\hat{\mathbf{V}}}_{\text{s}}^{\text{T}}\\ {\hat{\mathbf{V}}}_{\text{n}}^{\text{T}}\end{array}\right]\\ & ={\hat{\mathbf{V}}}_{\text{s}}{\hat{\mathbf{\Lambda }}}_{\text{s}}{\hat{\mathbf{V}}}_{\text{s}}^{\text{T}}+{\hat{\mathbf{V}}}_{\text{n}}{\hat{\mathbf{\Lambda }}}_{\text{n}}{\hat{\mathbf{V}}}_{\text{n}}^{\text{T}}.\end{array}$$ 由于信号是异步传输的，相当于每个用户的信号给相关矩阵${\hat{\mathbf{R}}}_{\text{z}}$贡献两个大特征值，因此${\hat{\mathbf{\Lambda }}}_{\text{s}}$的对角线元素包含${\hat{\mathbf{R}}}_{\text{z}}$的前$2K$个大特征值，${\hat{\mathbf{V}}}_{\text{s}}$的列向量张成信号子空间，记为$\text{Span}({\hat{\mathbf{V}}}_{\text{s}})$。
  第$k$个用户的${\mathbf{u}}_k^{\text{r}}$和${\mathbf{u}}_k^{\text{l}}$都位于信号子空间$\text{Span}({\hat{\mathbf{V}}}_{\text{s}})$中，那么${\mathbf{u}}_k={\mathbf{u}}_k^{\text{r}}+{\mathbf{u}}_k^{\text{l}}$也位于$\text{Span}({\hat{\mathbf{V}}}_{\text{s}})$中。${\mathbf{u}}_k$在$\text{Span}({\hat{\mathbf{V}}}_{\text{s}})$上的投影表示为 $$f({\mathbf{u}}_k)=({\mathbf{u}}_k^{\text{T}}{\hat{\mathbf{V}}}_{\text{s}})({\mathbf{u}}_k^{\text{T}}{\hat{\mathbf{V}}}_{\text{s}}{)}^{\text{T}}=\Vert {\mathbf{u}}_k^{\text{T}}{\hat{\mathbf{V}}}_{\text{s}}{\Vert }^2.$$ 第$k$个用户的延时的采样点数量$q_k$可以由下式获得 q ^ k = arg ⁡ max ⁡ q k ∈ { 0 , ⋯   , L N s − 1 } f ( u k ) . \hat{q}_k = \mathop{\arg\max}\limits_{q_k \in \{0,\cdots,LN_\text{s}-1\}}f(\boldsymbol{u}_k). q^​k​=qk​∈{0,⋯,LNs​−1}argmax​f(uk​). 上述最小化问题可以通过遍历$q_k$的所有可能值来找到最优值${\hat{q}}_k$来解决。 对于每个用户，上述方法需要$L{N}_{\text{s}}$次搜索。第$k$个用户的延时为${\hat{\tau }}_k={\hat{q}}_k{T}_{\text{s}}$。

function [delay,sigma] = subspace_geo_channel_estimation(rec_data,ss_code,sf,M,L_bit,K,b,sps,isShape)
% Subspace-based channel estimation, geometric approach
% rec_data 接收信号
% ss_code 扩频码
% sf 扩频因子
% M 快拍数
% L_bit 信息比特数
% K 用户数
% b 成型滤波抽头系数
% sps 上采样率
% isShape 发射信号是否成型滤波 1 滤波, 0 不滤波
if M >= L_bit
    M = L_bit-1;
end
y = zeros(sf*sps,M);
for m = 1:M
    for l = 1:sf*sps
        y(l,m) = rec_data(l+sf*sps*(m-1));
    end
end
Ry = (y*y')./M;
[V,D] = eig(Ry);
[D,ind] = sort(diag(D),'descend'); % 特征值按由大到小排
V = V(:,ind);
V_S = V(:,1:2*K); % signal subspace
if isShape == 1 % 成型滤波
    ss_code_upsample = upfirdn(ss_code',b,sps)';
    ss_code_upsample = ss_code_upsample(:,round(length(b)/2)-round(sps/2)+1:round(length(b)/2)-round(sps/2)+sps*sf);
else
    ss_code_upsample = rectpulse(ss_code',sps)';% 矩形成型
end
J = zeros(K,sf*sps);
for k = 1:K  
    for v = 0:sf*sps-1
        a_r_0 = [ss_code_upsample(k,sf*sps+1-v:end),zeros(1,sf*sps-v)]';
        a_l_0 = [zeros(1,v),ss_code_upsample(k,1:sf*sps-v)]';
        a_r = a_r_0;
        a_l = a_l_0;
        J(k,v+1) = norm((a_r+a_l)'*V_S)^2;
    end
end
[~,ind] = max(J,[],2);
v = ind-1;
delay = mod(v,sf*sps);
sigma = mean(D(2*K+1:end)); % 噪声方差
end

3、算法仿真

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

%% 仿真参数
date = '5_28';
if(~exist(['.\',date,'sim_dadta'],'dir'))
    mkdir(['.\',date,'sim_data']);
end
K = 3; % 用户数
P = 10 ; % 上采样率 samples/chip
isShape = 1; %是否成型滤波 1 滤波, 0 不滤波
isDC = 0; % 接收机直流耦合 1 直流耦合, 0 交流耦合, 可以直流耦合接收，后面在代码里去直
Target_User = 1; % 目标用户
noise_power = 22:-2:8; % noise power dBW
target_user_power = 0; % AC power (variance) dBW
interference_user_power = [0,5,10,20]; % AC power (variance) dBW
Times = 20;
win_size = [20,25,30,35,40,50:25: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
% 用户发送数据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)];    
L_ss = 2^(length(ss_polynomial)-1)-1; % length of spread spectrum pn sequence, spreading factor
L_bits = 300;
delay_array = 0:8:L_ss*P-1;
%% 生成发送数据
ss_code = zeros(K,L_ss);
user_bits = zeros(K,L_bits);
user_ss_data = zeros(K,L_bits*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_bits_upsample = rectpulse(user_bits(k,:),L_ss);
    user_ss_data(k,:) = user_bits_upsample.*repmat(ss_code(k,:),1,L_bits); 
end
% 上采样，成型滤波
if isShape == 1
    sps = P; % 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 = P;
    user_ss_data = rectpulse(user_ss_data',sps)';% 矩形成型
end
if isDC == 1 % 光信号为正的实信号
    user_ss_data = user_ss_data-min(user_ss_data,[],2);
end
user_ss_data = user_ss_data./vecnorm(user_ss_data,2,2); %能量归一化
user_ss_data = user_ss_data.*sqrt(length(user_ss_data));%功率归一化
clear user_bits_temp;
clear user_bits_upsample;
%% 接收
user_ss_data(Target_User,:) = user_ss_data(Target_User,:)*sqrt(10^(target_user_power/10));
L_data = length(user_ss_data);
L_sample = L_data+P*L_ss;
% 每行先固定一个snapshots遍历interference_user_power，再遍历snapshots
% 成功捕获时才计算rmse
acquisition_pro_corr = zeros(length(win_size)*length(interference_user_power),length(noise_power));
acquisition_rmse_corr = zeros(length(win_size)*length(interference_user_power),length(noise_power));
acquisition_pro_geo = zeros(length(win_size)*length(interference_user_power),length(noise_power));
acquisition_rmse_geo = zeros(length(win_size)*length(interference_user_power),length(noise_power));
for t = 1:Times
    for n = 1:length(noise_power)
        for d = 1:length(delay_array)
            if isShape == 1 % 延时参考
                delay_ref = round(length(b)/2)-round(sps/2)+[delay_array(d),0,0];
                delay_ref = mod(delay_ref,sps*L_ss);
            else
                delay_ref = [delay_array(d),0,0];
            end         
            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)
                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
                background_noise = wgn(1,L_sample,noise_power(n));% background noise
                rec_data = sum(send_data,1)+background_noise;
                if isDC == 1
                    rec_data = rec_data-mean(rec_data); % DC block
                end
                for m = 1:length(win_size)
                    M = win_size(m);
                    % sliding correlator
                    delay = corr_channel_estimation(rec_data,ss_code,L_ss,M,L_bits,K,b,sps,isShape);
                    if abs(delay(Target_User)-delay_ref(Target_User)) < P/2
                        acquisition_pro_corr((m-1)*length(interference_user_power)+p,n) = acquisition_pro_corr((m-1)*length(interference_user_power)+p,n)+1;
                        acquisition_rmse_corr((m-1)*length(interference_user_power)+p,n) = acquisition_rmse_corr((m-1)*length(interference_user_power)+p,n)+(delay(Target_User)-delay_ref(Target_User))^2;
                    end
                    % subspace-based geometric approach
                    [delay,~] = subspace_geo_channel_estimation(rec_data,ss_code,L_ss,M,L_bits,K,b,P,isShape);
                    if abs(delay(Target_User)-delay_ref(Target_User)) < P/2
                        acquisition_pro_geo((m-1)*length(interference_user_power)+p,n) = acquisition_pro_geo((m-1)*length(interference_user_power)+p,n)+1;
                        acquisition_rmse_geo((m-1)*length(interference_user_power)+p,n) = acquisition_rmse_geo((m-1)*length(interference_user_power)+p,n)+(delay(Target_User)-delay_ref(Target_User))^2;
                    end
                end
            end
        end
    end
end
acquisition_rmse_corr = sqrt(acquisition_rmse_corr./acquisition_pro_corr)./P;
acquisition_pro_corr = acquisition_pro_corr./(t*d);
acquisition_rmse_geo = sqrt(acquisition_rmse_geo1./acquisition_pro_geo1)./P;
acquisition_pro_geo = acquisition_pro_geo1./(t*d);
SNR = target_user_power-noise_power;
EbN0 = SNR-10*log10(1/L_ss)+10*log10(P/2);
MAI = interference_user_power-target_user_power; % near far ratio
save(['.\',date,'sim_data\',date,'sim_estimation_pro_rmse_avr'],'acquisition_pro_corr','acquisition_rmse_corr',...
      'acquisition_pro_geo','acquisition_rmse_geo','SNR','EbN0','MAI','win_size');

function p=PnCode(polynomial,reg)
%  PN码产生器函数
%　polynomial的长度＝reg的长度+1,polynomial的值不能为全0
%  polynomial为本原多项式，从左到右依次为高位到低位，且最高位与最低位必须为1;低位表示延时一个周期，高位依次顺延
%  reg为置寄存器初始值,也相当于PN码的初始相位，左边为高位，如[1 0 0 1 0]表示延时5个周期的寄器和2个周期的寄存器初值为1
%　如产生5级31位的PN码，则多项式形式为[1　* * * * 1]
%  例：5级PN码45E，参数为[1 0 0 1 0 1],左边为高位
%　PN：0 1 0 0 0 0 1 0 1 0 1 1 1 0 1 1 0 0 0 1 1 1 1 1 0 0 1 1 0 1 0
grade=length(polynomial)-1;%根据多项式计算延时级数
PN_Length=(2^grade-1);     %计算PN码一个周期的长度 
%找出本原多项式中除最低位外为1的位,并依次存放在寄存器c中
%例如对于ploynomial＝[1 0 0 1 0 1]，则c(1)=2,c(2)=5
n=0;                         
c=zeros(1,grade);
for i=grade:-1:1
    if polynomial(i)==1
        n=n+1;
        c(n)=grade+1-i;
    end
end  
p = zeros(1,PN_Length);
for i=1:PN_Length %从最高延时的寄存器中输出PN码
    p(i)=reg(1);
    m = mod(reg*polynomial(1:grade)',2);%完成各抽头寄存器取值的模2加
    reg(1:(grade-1)) = reg(2:grade);%寄存器的值依次移位
    reg(grade)=m;
end

  代码有点多，可能有的函数没贴上来，缺代码的话请留言、私信或者点此下载未删减的全套代码。
  捕获概率定义为估计的延时与实际延时的误差小于半个码片周期的概率 $$P\left(\left|{\tau }_k-{\hat{\tau }}_k\right|<\frac{T_{\text{c}}}{2}\right),$$ 其中，${\tau }_k$表示实际延时，${\hat{\tau }}_k$表示估计的延时。这里给一个信噪比为$-4$ dB，用户2与用户1的功率比为$20$ dB时，用户1的捕获概率随快拍数变化的仿真结果，如下图所示。

  滑动相关法的捕获概率一直很差。这是因为受MAI和远近效应的影响，滑动相关法极有可能将互相关峰出现的位置作为延时位置。基于子空间的延时估计算法能克服MAI带来的负面影响，它的捕获概率随快拍数的增加而增大，并且在快拍数为$75$时达到$100\%$的捕获概率。

参考文献

[1] PICKHOLTZ R, SCHILLING D, MILSTEIN L. Theory of spread-spectrum communications-a tutorial[J]. IEEE Transactions on Communications, 1982, 30(5): 855-884.

[2] BENSLEY S E, AAZHANG B. Subspace-based channel estimation for code division multiple access communication systems[J]. IEEE Transactions on Communications, 1996, 44(8):1009-1020.

[3] STROM E G, PARKVALL S, MILLER S L, et al. Propagation delay estimation in asynchronous direct-sequence code-division multiple access systems[J]. IEEE Transactions on Communications, 1996, 44(1):84-93.