基于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');  %在坐标區上添加圖例