BOC调制同步系统的鉴相器性能MATLAB仿真

UP目录

一、理论基础

二、核心程序

三、测试结果

---

一、理论基础

       随着全球导航卫星系统的高速发展,导航系统的数量也越来越多,比如使用最广泛的GPS导航系统,以及越来越备受关注的中国北斗导航系统等.因此导航频段变得越来越拥挤,且各个频段内的信号相互干扰,在如此情况下,一种二进制偏移载波(Binary Offset Carrier,BOC)调制信号被成功提出,用来解决该问题.而伴随着对BOC信号的深入研究,一系列BOC衍生信号也被提出来应用于各个导航系统中.此外,实际环境中存在着许多干扰因素,给BOC及其衍生信号的捕获带来了困难,因此也成为了国内外学者研究的热点与难点.本文在详细研究了BOC及其衍生信号结构和特性的基础上,着重研究了低信噪比下BOC信号,高动态环境下BOC和高阶双二进制偏移载波(Double Binary Offset Carrier,DBOC)调制信号的捕获问题,以及使用改进的方法实现BOC及其衍生信号的通用无模糊捕获.

        BOC调制是Jhon W.Betz最先提出的，他是Galileo信号设计的一个领军人物。结合Galileo的频谱设计，以L1频段为例，我们可以看到，由于国际电联给导航频段分配的限制，Galileo和GPS必须公用一个带宽，最理想的中心频点已经被GPS占了，并且GPS的C/A信号（BPSK(1)）也把中心频段占了，因此Galileo信号只能避开C/A信号，而把功率谱分裂成两个，就很好的解决了这个问题。我个人理解，BOC信号的设计实际是信号功率谱赋形的一种反推过程，是根据实际需要而生成的一种调制方式。
        BOC 调制就是在原有 BPSK 调制的基础上，再加上一个二进制副载波对 BPSK 信号进行二次扩频。由于方波易于生产，因此采用方波来替代正弦波，这样可以节约硬件资源。一般用 BOC（sf ,cf）来表示sf 代表副载波频率，cf 代表伪码速率。因为sf 和cf都是 1.023MHz 的整数倍，所以文献中用 BOC(m,n)的表示形式，其中 m 表示的是副载波频率，n 表示的是扩频码速率，它们分别表示 1.023MHz 的 m 倍和 n 倍。
       BOC 调制的频谱可以分为两个对称独立的边带，它的自相关函数是一种多峰结构，其中主峰较窄，所以有可能对接收机处理方式做一些改变，以便在实现时获得性能的改善。另外一种可行的方法是使鉴相器工作在一个带宽信号的相对函数上，进而保持对主峰的粗跟踪，之后再转入用双边带信号工作的鉴相器，来作精跟踪。由此可见，BOC 调制信号的频谱中上边带和下边带的冗余，自相关函数中的多峰，如果加以开发利用有可能在接收机的信号截获、码跟踪以及数据解调等处理中带来利益。

        原理是在原有的BPSK调制基础上，再增加一个二进制副载波（目前主要是由正弦或余弦型符号函数构成的副载波，即形似sgn(sin(t))或sgn(cos(t))，以正弦或余弦信号为参数的符号函数）。这种调制方式的最大特点是，其功率谱的主瓣分裂成对称的两部分，而且根据选择的参数不同，两个分裂主瓣的距离也可以变化。一般常用的表示方式为BOC(m,n)的形式，其中m表示的是副载波频率，n表示的是扩频码速率，具体数值分别是1.023MHz的m倍和n倍。

       关于BOC调制，Betz的文章《Binary Offset Carrier Modulations for Radionavigation》是比较经典的，特别是附录里关于功率谱表达式的推导，让我省却了很多自行推导的烦恼。不过里面有一处印刷错误，有兴趣的朋友应该不难发现。中文资料里我认为邱致和的《GPS M码信号的BOC调制》，可以作为参考，特别是关于BOC基本概念的介绍。

       MBOC－Multiplexed Binary Offset Carrier，从英文名可以看出这就是BOC副载波调制信号的一种复用方式。这是由Guenter W.Hein领导的GPS信号设计团队和Jhon W.Betz领导的Galileo信号设计团队共同提出的一种调制方式。目前经过优选，主要讨论和设计应用的是BOC(1,1)和BOC(6,1)的组合。具体根据数据通道和导频通道的功率分配要求，以及采取具体的调制方式不同，可以有多种组合，具体可参考相关文献。目前基本都是从BOC(1,1)和BOC(6,1)的功率分配角度来讨论的。

       MBOC只是一种信号复用的统称，其具体实现目前主要有两种，即CBOC（Composite BOC）和TMBOC（time-multiplexed BOC）。CBOC简单的说，是根据BOC(1,1)和BOC(6,1)不同的功率（幅值）权重构成的4电平符号来实现的调制，是幅值的复合式实现。而TMBOC则是一种类似时分复用的方式，即规定一组码片的长度，在这组码片里固定的几个位置里是BOC(6,1)，其他位置都是BOC(1,1)。两种方式都能满足功率谱分配的要求，但在功率谱谱形上还是有所不同的。

二、核心程序

..........................................................
Tcc           = [-1:0.1:1];
Per           = zeros(size(Tcc));

MTKL          = 1000;%蒙特卡洛循环次数
for i = 1:length(Tcc)
    error=0;
    for jj = 1:MTKL
    CN0			  =	40;
    LN            = 40;
    Scale         = 2048;
    f0            = 1.023e6;%基准频率
    BOCm          = 10;
    BOCn          = 5;        
    BW		      =	2*(BOCm+BOCn)*f0; 

    SNR			  =	CN0 - 10*log10(BW);	
    %符号速率为50sps
    Rb            = 50; 
    navdatLength  = 1;
    %仿真持续时间为10个数据比特
    D             = navdatLength/Rb; 
 
    %副载波频率
    fs            = 50*f0;
    %码速率
    fc            = 5*f0; 
    fIF           = 20*f0;
    fsamp         = 100e6;%采样频率

    Code_cc       = fc*2^LN/fsamp;
    SubCarrier_cc = fs*2^LN/fsamp;
    Carrier_cc    = fIF*2^LN/fsamp;

    %信号源
    Cd_accumulator  = 0;
    Sc_accumulator  = 0;
    S_accumulator   = 0;
    ROM_sin         = sin(2*pi/Scale*(0:Scale-1));
    ROM_code        = func_gold_gen(Scale/2-1);

    La_code         = 0;
    La_subcarrier   = 0;
    La_fll          = 0;
    La_pll          = 0;
    Lfcode          = Code_cc;
    Lfsubcarrier    = SubCarrier_cc;
    Lf_fll          = Carrier_cc;
    Lf_pll          = Carrier_cc;

    Nl_code         = 0;
    Nl_subcarrier   = 0;
    Nl_fll          = 0;
    Nl_pll          = 0;
    
    %相干时间,2ms
    Tcoh            = 2e-3; 
    %仿真跟踪时间
    T_track         = 4;     
    Lens            = round(Tcoh*fsamp);
    CAIndex         = floor(fsamp/fc);
    %码环的自然频率
    Nfreq           = 0.94;      
    %环路滤波器系数
    C1              = (sqrt(2)*Nfreq + Nfreq^2*Tcoh);
    C2              = (sqrt(2)*Nfreq);
    loopnum         = round(T_track/Tcoh);
    [signal_BOC]    = func_BOC_signal(2,Code_cc,SubCarrier_cc,Carrier_cc);

    for index = 1:400
        i
        jj
        index
        %信号跟踪
        [signal_BOC] = func_BOC_signal(Lens,Code_cc,SubCarrier_cc,Carrier_cc);
        %加高斯白噪声
        signal_BOC0  = awgn(signal_BOC,SNR,'measured');
        Lenss        = round(abs(200*Tcc(i)));
        if Tcc(i)>0
           signal_BOC1 = [signal_BOC0(Lenss:end),signal_BOC0(1:Lenss-1)];
        end
        if Tcc(i)==0
           signal_BOC1 = signal_BOC0;
        end        
        if Tcc(i)<0
           signal_BOC1 = [signal_BOC0(end-Lenss+1:end),signal_BOC0(1:end-Lenss)];
        end          
        
        
        %本地码载波产生
        [Lprn_E,Lprn_L,Lboc_E,Lboc_L,Fcos,Fsin,Pcos,Psin] = func_NCOprn(Lens);
        %前端滤波
        B          = fir1(64,0.75,'low');
        signal_BOC = filtfilt(B,1,signal_BOC1);
        
        %相关运算
................................................................

        %码鉴别滤波，计算频率字
        pow_E                  = Dcm_I_E^2+Dcm_Q_E^2;
        pow_L                  = Dcm_I_L^2+Dcm_Q_L^2;
        code_dis               = (pow_E-pow_L)/(pow_E+pow_L)/2;
        code_dis_result(index) = code_dis;

        if index ==1
            last_out_DDL = 0;
            last_in_DDL  = 0;
        end
        code_filter               = last_out_DDL + [code_dis last_in_DDL]*[C1 (-1)*C2]';
        last_out_DDL              = code_filter;
        last_in_DDL               = code_dis;
        code_filter_result(index) = code_filter;
         
        Lfcode          = (fc+code_filter)*2^LN/fsamp;
        
         
        Lfsubcarrier    = (fs+code_filter)*2^LN/fsamp;
    end
    error(jj) = mean(code_filter_result(1:end));
    end
    Per(i) = mean(error);
end

figure;
plot(Tcc,Per,'b-o');
grid on
xlabel('input Code Delay (chips)');
ylabel('Discriminator output');
grid on
up74

三、测试结果

通过matlab2021a仿真结果如下：