基于BOC的信号捕获

1.问题描述：

 随着全球导航卫星系统的高速发展,导航系统的数量也越来越多,比如使用最广泛的GPS导航系统,以及越来越备受关注的中国北斗导航系统等.因此导航频段变得越来越拥挤,且各个频段内的信号相互干扰,在如此情况下,一种二进制偏移载波(Binary Offset Carrier,BOC)调制信号被成功提出,用来解决该问题.而伴随着对BOC信号的深入研究,一系列BOC衍生信号也被提出来应用于各个导航系统中.此外,实际环境中存在着许多干扰因素,给BOC及其衍生信号的捕获带来了困难,因此也成为了国内外学者研究的热点与难点.本文在详细研究了BOC及其衍生信号结构和特性的基础上,着重研究了低信噪比下BOC信号,高动态环境下BOC和高阶双二进制偏移载波(Double Binary Offset Carrier,DBOC)调制信号的捕获问题,以及使用改进的方法实现BOC及其衍生信号的通用无模糊捕获.论文的主要工作如下:(1)针对低信噪比条件下BOC信号的捕获问题,研究了一种基于批处理和Teager-Kaiser(TK)算子的联合捕获算法.该方法首先对接收信号和组合扩频码进行批处理和平均处理,并且将接收信号转换到频域进行多次循环移位,达到多普勒频偏补偿的效果;然后将多个频偏补偿序列与组合扩频码进行圆周相关运算;最后将其中的最大相关值序列经过TK算子处理.研究表明,在同一条件下,本文方法有效提升了检测性能,

2.部分程序：

 

longSignal = data;
%% Initialization =========================================================
% Find number of samples per spreading boc code period
samplesPerCode = round(settings.samplingFreq * (settings.BOCcodeLength/ ...
                        settings.BOCcodeFreq));     %每个伪码周期的采样点数

% Create two 1msec vectors of data to correlate with and one with zero DC
signal1 = longSignal(1 : samplesPerCode);           %第一个伪码周期对应的采样数据
signal2 = longSignal(samplesPerCode+1 : 2*samplesPerCode);  %第二个伪码周期对应的采样数据

signal0DC = longSignal - mean(longSignal);         %去直流

% Find sampling period
ts = 1 / settings.samplingFreq;                  %采样周期

% Find phase points of the local carrier wave in one spreading code period
phasePoints = (0 : (samplesPerCode-1)) * 2 * pi * ts;     %对应的采样相位

% Number of the frequency bins for the given acquisition band (500Hz steps)
numberOfFrqBins = round(settings.acqSearchBand * 2) + 1;   %每次搜索步长为500hz，总搜索次数

% Generate all C/A & BOC codes and sample them according to the sampling freq.
[caCodesTable,bocCodesTable] = makeCaTable(settings); %%改caCodesTable = makeCaTable(settings);
                                                       %使用makeCaTable函数产生所有
                                                       %的BOC序列矩阵与C/A码矩阵
%--- Initialize arrays to speed up the code -------------------------------
% Search results of all frequency bins and code shifts (for one satellite)
results     = zeros(numberOfFrqBins, samplesPerCode);   %结果为多普勒频率和码序列二维搜索过程

% Carrier frequencies of the frequency bins
frqBins     = zeros(1, numberOfFrqBins);               %载波频率初始化

%--- Initialize acqResults ------------------------------------------------
% Carrier frequencies of detected signals
acqResults.carrFreq     = zeros(1, 32);
% C/A code phases of detected signals
acqResults.codePhase    = zeros(1, 32);
% Correlation peak ratios of the detected signals
acqResults.peakMetric   = zeros(1, 32);               %初始化捕获结果

fprintf('(');
tic;
%%

% Perform search for all listed PRN numbers ...
%for PRN = settings.acqSatelliteList                %图6.8的并行码相位空间捕获过程（双通道比较？）

%% Correlate signals ======================================================   
    %--- Perform DFT of C/A code ------------------------------------------
     PRN=1;
    bocCodeFreqDom1 = conj(fft(bocCodesTable(PRN, :))); %%改caCodeFreqDom = conj(fft(caCodesTable(PRN, :)));
    bocCodeFreqDom2 = conj(fft(caCodesTable(PRN, :)));                                                   %本地接收码
                                                       %PRN码发生器经傅里叶变换再经共轭
   
%     Raspect=(abs(Rboc)).^2-(abs(Rboc_ca)).^2 ; 
    %--- Make the correlation for whole frequency band (for all freq. bins)
    for frqBinIndex = 1:numberOfFrqBins

        %--- Generate carrier wave frequency grid (0.5kHz step) -----------
        frqBins(frqBinIndex) = settings.IF - ...
                               (settings.acqSearchBand/2) * 1000 + ...
                               0.5e3 * (frqBinIndex - 1);   %不同多普勒频移下的载波频率

        %--- Generate local sine and cosine -------------------------------
        sinCarr = sin(frqBins(frqBinIndex) * phasePoints);
        cosCarr = cos(frqBins(frqBinIndex) * phasePoints);   %本地振荡器

        %--- "Remove carrier" from the signal -----------------------------
        I1      = sinCarr .* signal1;
        Q1      = cosCarr .* signal1;
        I2      = sinCarr .* signal2;
        Q2      = cosCarr .* signal2;                 %去载波

        %--- Convert the baseband signal to frequency domain --------------
        IQfreqDom1 = fft(I1 + j*Q1);
        IQfreqDom2 = fft(I2 + j*Q2);                  %fft

        %--- Multiplication in the frequency domain (correlation in time
        %domain)
        convCodeIQ11 = IQfreqDom1 .* bocCodeFreqDom1;%%改caCodeFreqDom;
        convCodeIQ12 = IQfreqDom1 .* bocCodeFreqDom2;
%         convCodeIQ1=(abs( convCodeIQ11)).^2-(abs( convCodeIQ12)).^2 ;
        convCodeIQ21 = IQfreqDom2 .* bocCodeFreqDom1;%%caCodeFreqDom;
        convCodeIQ22 = IQfreqDom2 .* bocCodeFreqDom2;
%         convCodeIQ2=(abs( convCodeIQ21)).^2-(abs( convCodeIQ22)).^2 ;
        %--- Perform inverse DFT and store correlation results ------------
        acqRes1 = abs(ifft(convCodeIQ11)) .^ 2-abs(ifft(convCodeIQ12)) .^ 2;
        acqRes2 = abs(ifft(convCodeIQ21)) .^ 2-abs(ifft(convCodeIQ22)) .^ 2;         %逆fft取平方
        
        %--- Check which msec had the greater power and save that, will
        %"blend" 1st and 2nd msec but will correct data bit issues
        if (max(acqRes1) > max(acqRes2))
            results(frqBinIndex, :) = acqRes1;
        else
            results(frqBinIndex, :) = acqRes2;
        end                                            %1.比较大小经过循环可以找到多普勒频率
        % results(frqBinIndex, :) = acqRes1;              %以及含最大的自相关的输出结果
                                                       %2.此时只循环了一个PRN码即只对一个可见
                                                       %星的不同频率进行了搜索
                                                       %3.生成每个频率所对应的结果矩阵                                              
    end % frqBinIndex = 1:numberOfFrqBins
    toc;
    figure(18);
    clf(18)
    phase=0 : (samplesPerCode-1);
    [X,Y]=meshgrid(phase*settings.BOCcodeLength/samplesPerCode,frqBins);
    surf(X,Y, results);

    title (['第(' int2str(PRN) ')颗可见星的ASPeCT算法捕获结果']);
    ylabel('载波频率 (MHz)'); xlabel('相位（码片）');zlabel('幅度');
        axis tight;
        
  
                                                       %***捕获效果图*********
        %会产生有两个峰的原因？                                             
%% Look for correlation peaks in the results ==============================
    % Find the highest peak and compare it to the second highest peak
    % The second peak is chosen not closer than 1 chip to the highest peak
    
    %--- Find the correlation peak and the carrier frequency --------------
    [peakSize, frequencyBinIndex] = max(max(results, [], 2)); %找到矩阵最大值所在行数

    %--- Find code phase of the same correlation peak ---------------------
    [peakSize codePhase] = max(max(results));         %找到矩阵最大值以及最大值的列数

     figure(21);
    clf(21)
    plot((codePhase-370:codePhase+370)*settings.BOCcodeLength/samplesPerCode, results(frequencyBinIndex, codePhase-370:codePhase+370))
    title (['第(' int2str(PRN) ')颗可见星的ASPeCT算法捕获结果']);
    xlabel('相位（码片）');ylabel('幅度');
    %--- Find 1 chip wide BOC code phase exclude range around the peak ----
    %samplesPerBOCCodeChip   = round(settings.samplingFreq / settings.BOCcodeFreq);%%%BOC code freq: the code after CA code modulating subcarrier
    samplesPerCACodeChip   = round(settings.samplingFreq / settings.codeFreq);%%CA code freq
    excludeRangeIndex1 = codePhase - samplesPerCACodeChip;%如保护带为左右各一个BOC码片（samplesPerBOCCodeChip），则可能导致第二个互相关峰为主/副峰，使捕获失败
    excludeRangeIndex2 = codePhase + samplesPerCACodeChip;%设定相关函数的码片延迟范围为以最大值为中心左右各一个码片

    %--- Correct C/A code phase exclude range if the range includes array
    %boundaries%[codePhase:samplesPerCode,1:codePhase-1]范围去掉第一峰的保护带（左右各一个CA
    %码片，即+/-samplesPerCACodeChip）得到的范围用以搜索第二个互相关峰。这时涉及到“跨界”问题
    if excludeRangeIndex1 < 1
        codePhaseRange = excludeRangeIndex2 : ...
                         (samplesPerCode + excludeRangeIndex1);   %上边带左限小于1时，转换到下一个周期
                         
    elseif excludeRangeIndex2 >= samplesPerCode
        codePhaseRange = (excludeRangeIndex2 - samplesPerCode) : ...
                         excludeRangeIndex1;                      %下边带右限超出范围时，转换到前一个周期
    else
        codePhaseRange = [1:excludeRangeIndex1, ...
                          excludeRangeIndex2 : samplesPerCode];
    end                                                           % codePhaseRange即为去掉最大峰值左右一个码片的数据以便找到第二峰值点

    %--- Find the second highest correlation peak in the same freq. bin ---
    secondPeakSize = max(results(frequencyBinIndex, codePhaseRange));

    %--- Store result -----------------------------------------------------
    acqResults.peakMetric(PRN) = peakSize/secondPeakSize;         %获取峰值比
    
    % If the result is above threshold, then there is a signal ...
    if (peakSize/secondPeakSize) > settings.acqThreshold 

%% Fine resolution frequency search =======================================
        
        %--- Indicate PRN number of the detected signal -------------------
        fprintf('%02d ', PRN);                                   %峰值比大于阈值，输出捕获的PRN号码，阈值为峰值比？
        
        %--- Generate 10*CA periods long C/A codes sequence for given PRN --------
        caCode1023 = generateCAcode(PRN);        
        codeValueIndex = floor((ts * (1:10*samplesPerCode)) / ...
                               (1/settings.codeFreq));
                           
        longCaCode = caCode1023((rem(codeValueIndex, settings.CAcodeLength) + 1));
        clear caCode1023
        %--- Generate subcarrier ----------------------------- 
        subCarr = square(settings.subcarrFreq * ts * (1:10*samplesPerCode) * 2 * pi);%sign(sin(settings.subcarrFreq * phasePoints));
%         subCarr = square(settings.subcarrFreq * ts *(codePhase:(codePhase + 10*samplesPerCode-1)) * 2 * pi);
        longBocCode = longCaCode.*subCarr;                        %产生该可见星本地的BOC码数据用于解调
        clear longCaCode subCarr
%         caCode1023 = generateCAcode(PRN);
%         [caCode, bocCode] = CodeExtend(caCode1023,settings.BOCm,settings.BOCn);
%         
%         codeValueIndex = floor((ts * (1:10*samplesPerCode)) / ...
%                                (1/settings.BOCcodeFreq));
%                            
%         longBocCode = bocCode((rem(codeValueIndex, settings.BOCcodeLength) + 1));
        
        %--- Remove BOC code modulation from the original signal ----------
        % (Using detected C/A code phase)
        xCarrier = ...
            signal0DC(codePhase:(codePhase + 10*samplesPerCode-1)) ...
            .* longBocCode;%longCaCode;                              %去BOC调制得到原始数据
        
        %--- Find the next highest power of two and increase by 8x --------
        fftNumPts = 8*(2^(nextpow2(length(xCarrier))));             %fft的点数
        
        %--- Compute the magnitude of the FFT, find maximum and the
        %associated carrier frequency 
        fftxc = abs(fft(xCarrier, fftNumPts));                     %做fft
        
        uniqFftPts = ceil((fftNumPts + 1) / 2);
        [fftMax, fftMaxIndex] = max(fftxc(5 : uniqFftPts-5));        %找最大值序号，其中5？
        
        fftFreqBins = (0 : uniqFftPts-1) * settings.samplingFreq/fftNumPts;
        
        %--- Save properties of the detected satellite signal -------------
        acqResults.carrFreq(PRN)  = fftFreqBins(fftMaxIndex);       
        acqResults.codePhase(PRN) = codePhase;                        %输出捕获结果
    
    else
        %--- No signal with this PRN --------------------------------------
        fprintf('. ');
    end   % if (peakSize/secondPeakSize) > settings.acqThreshold
%     end    % for PRN = satelliteList

%=== Acquisition is over ==================================================
fprintf(')\n');
 

3.仿真结论：

D-48