基于MATLAB编写的GNSS_SDR（GNSS软件接收机）——自学笔记（2）

上一篇中，由于初始化中涉及到了probeData()函数，因此本文就先了解一下probeData()函数，然后是重点，介绍了捕获、捕获后的绘图函数。

目录

• probeData(varargin)函数
• acquisition.m 捕获函数
• • 捕获函数的原理与功能
  • Initialization 初始化
  • 捕获过程
• plotAcquisition.m

probeData(varargin)函数

函数功能：通过该函数，绘制原始数据信息:时域图，频域图和直方图。

该函数能通过两种方式进行调用probeData(settings)或者probeData(fileName, settings)。

其中输入参数为:

fileName - name of the data file. File name is read from settings if parameter fileName is not provided.数据文件的名称，如果没有提供参数fileName，则从设置中读取文件名。

settings - receiver settings. Type of data file, sampling frequency and the default filename are specified here. 接收机参数设置，这里指定了数据文件的类型、采样频率和默认文件名。

由于存在两种调用probeData(varargin)函数的方式，因此首先需要检查参数的数量，使用nargin来判断输入变量个数。

function probeData(varargin)
%% Check the number of arguments 检查参数的数量=======================
if (nargin == 1)  %nargin是用来判断输入变量个数的函数，这样就可以针对不同的情况执行不同的功能。
    settings = deal(varargin{1}); %deal函数的作用：将单个输入值据赋值给所有输出参数。
    fileNameStr = settings.fileName;
elseif (nargin == 2)
    [fileNameStr, settings] = deal(varargin{1:2}); %分别将两个输入参数赋给fileNameStr, settings
    if ~ischar(fileNameStr)
        error('File name must be a string');
    end
else
    error('Incorect number of arguments');
end
           

接下来，通过对文件的读操作，将fileNameStr文件读入data中，并根据原始数据画时域图，频域图和直方图，程序注释如下：

%% Generate plot of raw data ==============================================
[fid, message] = fopen(fileNameStr, 'rb');  %如果打开不成功，则返回与系统相关的错误消息

if (fid > 0)   %表明打开成功
    % Move the starting point of processing. Can be used to start the
    % signal processing at any point in the data record (e.g. for long
    % records).移动处理的起点，可在数据记录的任何点开始信号处理(例如长记录)。
    
    fseek(fid, settings.skipNumberOfBytes, 'bof');    
%     fseek函数用于定位文件位置指针，其调用格式为：
%     status=fseek(fid, offset, origin)
%     其中fid为文件句柄，offset表示位置指针相对移动的字节数，origin表示位置指针移动的参照位置
%       offset> 0    Move toward the end of the file.
%       offset= 0    Do not change position.
%       offset< 0    Move toward the beginning of the file.
%       'bof' or -1  Beginning of file
%       'cof' or 0   Current position in file
%       'eof' or 1   End of file  


% Find number of samples per spreading code 找出每1ms（C/A码一周期为1ms）的采样点数
    samplesPerCode = round(settings.samplingFreq / ...
                           (settings.codeFreqBasis / settings.codeLength));
                      %理解为：采样频率/（C/A码的码率[Hz]/码片长度）
                      %按初始化参数计算=38.192MHz/（1.023*10^6/1023）
                      %也就可以理解为每1ms的采样点数 参考《GPS原理与接收机设计》P260，数模转换
                      
    % Read 10ms of signal
    [data, count] = fread(fid, [1, 10*samplesPerCode], settings.dataType);
    % 该语句利用Fread函数，从指定文件中读取二进制数据并写入矩阵data，Count参数用于返回成功读入的元素数量，为可选参数。
    
%       1、fid 文件标识 
%       2、sizeA     输出数组的维度
%           有3种参数,Inf、n、[m,n]，[m,n] 代表按列向量排列的m行n列的矩阵,n可以取Inf，但m不可以
%       3、precision 需要读取数据的类型和大小，默认'uint8=>double'
%           常见有uint,uint8、uint16等数据格式，需要根据源数据来确定

    fclose(fid);
    
    if (count < 10*samplesPerCode)
        % The file is to short
        error('Could not read enough data from the data file.');
    end
    
    %--- Initialization ---------------------------------------------------
    figure(100);
    clf(100);
    
    timeScale = 0 : 1/settings.samplingFreq : 5e-3;    
    
    %--- Time domain plot -------------------------------------------------
    subplot(2, 2, 1);
    plot(1000 * timeScale(1:round(samplesPerCode/50)), ...
         data(1:round(samplesPerCode/50)));
     
    axis tight;
    grid on;
    title ('Time domain plot');
    xlabel('Time (ms)'); ylabel('Amplitude');
    
    %--- Frequency domain plot --------------------------------------------
    subplot(2,2,2);
    pwelch(data-mean(data), 16384, 1024, 2048, settings.samplingFreq/1e6)
    %pwelch（）函数利用Welch'平均功率图法返回信号x的功率谱密度（PSD）
    
    axis tight;
    grid on;
    title ('Frequency domain plot');
    xlabel('Frequency (MHz)'); ylabel('Magnitude');
    
    %--- Histogram 直方图 -------------------------------------------------
    subplot(2, 2, 3.5);
    hist(data, -128:128)
    
    dmax = max(abs(data)) + 1;
    axis tight;
    adata = axis;
    axis([-dmax dmax adata(3) adata(4)]);
    grid on;
    title ('Histogram'); 
    xlabel('Bin'); ylabel('Number in bin');
else
    %=== Error while opening the data file ================================
    error('Unable to read file %s: %s.', fileNameStr, message);
end % if (fid > 0)
           

运行该函数可以得到结果图，如图：

acquisition.m 捕获函数

该函数的功能是：对收集到的“数据”执行冷启动采集，并搜索所有卫星的GPS信号，在参数结构体中的“acqsatellite elist”字段中列出。同时返回一个“acqResults”结构体，其中保存了检测信号的码相位和频率。

function acqResults = acquisition(longSignal, settings)

Inputs:

longSignal - 前11ms原始信号

settings - 接收机参数设置，能够提供有关采样频率和中频及其他参数的资料，包括拟获取的卫星清单。

Outputs:

acqResults - 函数在“acqResults”结构中保存检测信号的码相位和频率。如果对给定的卫星PRN没有检测到信号，字段“carrFreq”被设置为0。

捕获函数的原理与功能

捕获函数的功能是：

捕获函数每隔0.5kHz搜索一次GPS信号，在每个顿率搜索阶段，同时搜索码相位。每次搜索之后保存相关结果并进入下一频率单元，由此函数遍历所有的频段(用户定义的多普勒空间)，接下来寻找函数最大相关值(相关峰)。之后，寻找同一频率单元中的第二大相关值。然后计算其比值并作为信号检测的规则，该比值用于和接收机中预设的变量值acq_threshold进行比较。

检测器与采样频率无关，因此也与峰值大小及噪声电平无关。

如果第一峰值与第二峰值的比值超过门]限值，则利用FFT方法找到精确的载波频率，以协助跟踪环中的PLL启动信号跟踪。

5kHz的频率精度对于PLL跟踪太粗略了，因此需要在找到信号后对信号利用DFT细化频率估计。 （加长相干积分时间和非相干积分时间可以提高信噪比，讨论过程间《GPS原理与接收机设计》P361）

接收机对二维搜索的基本操作：

它们两者都复制一定频率的载波和一定相位的C/A码，然后经过混频和相关来检测复制信号与接收信号之间的相关程度，信号捕获过程中的复制信号参数值是根据搜索策略设置的。只有当接收机内部所复制的载波和C/A码信号与接收信号相一致时，相关器的输出功率才会达到最大值，而信号的捕获正是通过检查相关器的输出功率在何种复制载波频率和C/A码相位的下达到最大来实现的。如果最大输出功率超过信号捕获门限值，那么相应于该最大功率的复制信号参数值就是信号捕获对当前接收信号的参数估计值。

该例程使用的信号搜索捕获算法是：并行码相位搜索法

通过之前对信号捕获的一些基本概念和总体框架的介绍，可知信号捕获对应于各个不同伪码、频率和码相位的三维搜索单元上进行信号搜索，而搜索的主要操作是进行C/A码相关运算。根据相关运算实现方式的不同，信号搜索捕获算法可主要分为线性搜索、并行频率搜索和并行码相位搜索三种。

并行码相位搜索法：如下图所示，当数字中频信号分别与I支路和Q支路上某一频率的复制正弦和复制余弦载波信号混频后，并行码相位搜索捕获算法并不是让这些混频结果i和q通过数字相关器直接与复制C/A码做相关运算，而是对复数形式的混频结果i+ jq进行傅里叶变换，然后将变换结果与复制C/A码傅里叶变换的共轭值相乘，接着将所得的乘积经傅里叶反变换得到在时域内的相关结果，最后对这些相关值进行检测来判断信号是否存在。我们稍后将会证明这个傅里叶反变换结果等于混频结果信号i+ jq与复制CIA码信号在当前搜索频带上各个码相位处的相关值。在完成了对当前频带的搜索与检测后，接收机接着让载波数控振荡器复制另一个频率值的正弦和余弦载波，然后类似地完成对下一个频带的搜索与检测。在对同一个卫星信号不同频带内的搜索过程中，复制CIA码的相位可保持不变，相应地其傅里叶变换及其共轭值也保持不变。当搜索另-个卫星信号时，接收机可让CIA码发生器复制相应的另一个CIA码，然后重复上述在各个不同频带中的信号搜索过程。

函数参数为初始的数据记录、预先产生的C/A码表和settings结构体。Settings结构中包含的与捕获有关的变量如下:

(1) Acq_ satelliteList:指定一 组卫星伪随机序列。捕获只对特定卫星进行。默认状态为空表，启动搜索所有可见卫星1~32。

(2) Acq_searchBand:指定卫星信号搜索的频率范围，为kHz的整数倍,并以IF为中心频率。用于捕获函数的步长为0.5kHz。

(3) Acq_threshold:信号检测器的门限。

输出是数组结构变量acqResult,包含Acq_ satelliteList中指定的所有卫星的搜索结果。如果某一颗卫星的信号被检测到,那么对应的signalDetected会被置1。

分段解读该代码：

Initialization 初始化

首先，该文件需要调用makeCaTable.m文件，是生成C/A码的函数，下一节将会对该函数做进一步阅读。

其次，该程序中的一些参数比较难懂，总结如下：

samplesPerCode 每1码周期内（C/A码一周期为1ms）的采样点

signal1 %取第一个码周期（1ms）内，原始信号的采样数据

signal0DC %经过去除直流信号处理的原始信号数据，有11ms

ts %采样时间

phasePoints %产生一个码片周期内的与原始信号同样多的相位点

numberOfFrqBins %f_unc/f_bin*2+1=频带数目

caCodesTable %本地生成的C/A码

caCodeFreqDom %首先由图13.10可知，需要对C/A码发生器进行傅里叶变换，然后求复数共轭，得到该矩阵

frqBins(frqBinIndex) %以0.5kHz为步长的频率估计值

codePhase %相关值最大时的相位值

samplesPerCodeChip %1码片宽的采样点数 值大概是38左右

acqResults.peakMetric %最大相关和次大相关的比值

acqResults.carrFreq %追踪载频的结果

初始化程序注释如下

function acqResults = acquisition(longSignal, settings)
%% Initialization =========================================================
% Find number of samples per spreading code
% 求解每1码片周期内（C/A码一周期为1ms）的采样点数
samplesPerCode = round(settings.samplingFreq / ...
                        (settings.codeFreqBasis / settings.codeLength));

% Create two 1msec vectors of data to correlate with and one with zero DC
% 创建两个1ms的相关数据向量和一个零直流数据
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 产生相位用的
phasePoints = (0 : (samplesPerCode-1)) * 2 * pi * ts;

% Number of the frequency bins for the given acquisition band (500Hz steps)
% 给定采集带宽的频带数量(500Hz步长)

%计算过程： 
% f_bin ： 对载波频率的搜索步长（又称，频带宽度）
% f_unc ： 信号频率的不定值  《设计》P351
% f_unc/f_bin = settings.acqSearchBand

numberOfFrqBins = round(settings.acqSearchBand * 2) + 1;%f_unc/f_bin*2+1=频带数目

% Generate all C/A codes and sample them according to the sampling freq.
% 生成所有的C/A代码，并根据采样频率进行采样。
caCodesTable = makeCaTable(settings);  %需要调用include文件夹下的makeCaTable.m文件


%--- 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);  %构建一个1行numberOfFrqBins列的矩阵，并赋初值


%--- Initialize acqResults 初始化捕获结果-----------------------------------
% Carrier frequencies of detected signals 检测信号的载波频率
acqResults.carrFreq     = zeros(1, 32);
% C/A code phases of detected signals  检测信号的C/A码相位
acqResults.codePhase    = zeros(1, 32);
% Correlation peak ratios of the detected signals  被检测信号的相关峰值比
acqResults.peakMetric   = zeros(1, 32);

fprintf('(');
           

捕获过程

此过程使用并行吗相位的搜索捕获算法，原理在上文已经阐述，下面为程序注解：

for PRN = settings.acqSatelliteList
%% Perform search for all listed PRN numbers ...执行搜索所有列出的PRN卫星
    %% Correlate signals 相关信号=========================================
    %--- Perform DFT of C/A code 计算CA码的相关-----------------------------
    caCodeFreqDom = conj(fft(caCodesTable(PRN, :)));  %首先由图13.10可知，需
                                                      %要对C/A码发生器进行傅
                                                      %里叶变换，然后求复数共轭
    %--- Make the correlation for whole frequency band (for all freq. bins)
    % 对整个频带进行相关性分析
    for frqBinIndex = 1:numberOfFrqBins

        %--- Generate carrier wave frequency grid (0.5kHz step) -----------
        % 等间距产生载波(0.5kHz步长)
        frqBins(frqBinIndex) = settings.IF - ...
                               (settings.acqSearchBand/2) * 1000 + ...
                               0.5e3 * (frqBinIndex - 1);
                           % （好像没问题了）？？？应该不能除以2吧？为啥要除以2？或者把
                           % （暂时保留）settings.acqSearchBand换成numberOfFrqBins也行

        %--- Generate local sine and cosine 产生了sin and cos---------------
        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);

        %--- Multiplication in the frequency domain (correlation in time
        %domain)  频域相乘(时域相关)
        convCodeIQ1 = IQfreqDom1 .* caCodeFreqDom;  %两者通过乘法器相乘
        convCodeIQ2 = IQfreqDom2 .* caCodeFreqDom;

        %--- Perform inverse DFT and store correlation results ------------
        % 执行逆DFT并存储相关结果，acqRes1，acqRes2均为矩阵形式
        acqRes1 = abs(ifft(convCodeIQ1)) .^ 2;
        acqRes2 = abs(ifft(convCodeIQ2)) .^ 2;
        
        %--- 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
        
    end % frqBinIndex = 1:numberOfFrqBins

%% 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
    % 找到最高的互相关峰值，并把它和第二高峰做比较，第二个峰值选择不小于1个码片的最高峰值
    %--- 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));

    %--- Find 1 chip wide C/A code phase exclude range around the peak ----
    % 要计算1码片以外的其他情况的互相关峰值
    % samplesPerCodeChip：1码片宽的采样点数 值大概是38左右
    samplesPerCodeChip   = round(settings.samplingFreq / settings.codeFreqBasis);
    excludeRangeIndex1 = codePhase - samplesPerCodeChip;
    excludeRangeIndex2 = codePhase + samplesPerCodeChip;

    %--- Correct C/A code phase exclude range if the range includes array
    %boundaries 如果取值在边界附近，需要修正C/A码相位
    if excludeRangeIndex1 < 2
        codePhaseRange = excludeRangeIndex2 : ...
                         (samplesPerCode + excludeRangeIndex1);
                         
    elseif excludeRangeIndex2 >= samplesPerCode
        codePhaseRange = (excludeRangeIndex2 - samplesPerCode) : ...
                         excludeRangeIndex1;
    else
        codePhaseRange = [1:excludeRangeIndex1, ...
                          excludeRangeIndex2 : samplesPerCode];
    end

    %--- 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 更好地解决频率搜索问题============
% 由于0.5kHz的频率精度对于PLL跟踪太粗略了，因此需要采用下面的方法更好地搜索频率
% 其原理是增大非相干采样的长度，从而增加精度。
        %--- Indicate PRN number of the detected signal -------------------
        fprintf('%02d ', PRN);
        
        %--- Generate 10msec long C/A codes sequence for given PRN --------
        caCode = generateCAcode(PRN);
        
        codeValueIndex = floor((ts * (1:10*samplesPerCode)) / ...
                               (1/settings.codeFreqBasis));  %floor取整函数
              
        longCaCode = caCode((rem(codeValueIndex, 1023) + 1));  %rem取余数
        %产生一个10ms的C/A码的矩阵，其采样频率已知
        %--- Remove C/A code modulation from the original signal ----------
        % (Using detected C/A code phase)
        xCarrier = ...
            signal0DC(codePhase:(codePhase + 10*samplesPerCode-1)) ...
            .* longCaCode;
        
        %--- Find the next highest power of two and increase by 8x --------
        fftNumPts = 8*(2^(nextpow2(length(xCarrier))));
        
        %--- Compute the magnitude of the FFT, find maximum and the
        %associated carrier frequency 
        fftxc = abs(fft(xCarrier, fftNumPts)); 
        
        uniqFftPts = ceil((fftNumPts + 1) / 2);
        [fftMax, fftMaxIndex] = max(fftxc(5 : uniqFftPts-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');
           

plotAcquisition.m

绘制捕获函数的结果图，程序比较简单，实现了函数绘制了条形结果图，且没有显示未在清单上的卫星图。结果图如图所示：

程序简要注释如下：

function plotAcquisition(acqResults)
%% Plot all results =======================================================
figure(101);

hAxes = newplot();  %newplot 为后续图形命令准备好图窗和坐标区。
                    %为后续图形命令准备图窗和坐标区并返回当前坐标区。

bar(hAxes, acqResults.peakMetric); % 将图形绘制到 hAxes 指定的坐标区中

title (hAxes, 'Acquisition results');
xlabel(hAxes, 'PRN number (no bar - SV is not in the acquisition list)');
ylabel(hAxes, 'Acquisition Metric');

oldAxis = axis(hAxes);
axis  (hAxes, [0, 33, 0, oldAxis(4)]);
set   (hAxes, 'XMinorTick', 'on');
set   (hAxes, 'YGrid', 'on');

%% Mark acquired signals标记捕获的信号 ===============================

acquiredSignals = acqResults.peakMetric .* (acqResults.carrFreq > 0);

hold(hAxes, 'on');
bar (hAxes, acquiredSignals, 'FaceColor', [0 0.8 0]);
hold(hAxes, 'off');

legend(hAxes, 'Not acquired signals', 'Acquired signals');  %在坐标区上添加图例