大学毕业设计一席谈之三十四 GPS卫星采样信号的跟踪(1) 展示程序和跟踪结果

知识的学习和能力提高过程是漫长的，也是循序渐进的！请看完上述付费系列文章，再来学习本系列文章，否则学起来太吃力！先来看看码跟踪环路的基础知识，然后再回顾锁相环的知识。有了这两个环路知识，才能对信号进行跟踪解调。

至上而下对图进行说明。深度为20的FIFO用于存储接收数据，由高速时钟驱动，每次从乒乓buffer读取4个4.092MHz采样频率的零中频数据，相当于1个码片时间长度的数据。数据通道对应FIFO的抽头，5路数据对应5个抽头，每个抽头映射到FIFO内的一个数据，中间抽头映射到FIFO第11+TapOffset个数据（见图），5个抽头之间的相互距离可以调整，间隔由控制字Tap_Width控制。比如在粗跟踪过程中，相互距离为2个采样数据，即1/2码片长；精跟踪过程中，相互距离为1个采样数据，即1/4码片长；惯性捕获时，相互距离为4个采样数据，即1个码片。本地码生成器主要产生0~1022个码片，方案采用查找表存储C/A码，由caIndex值寻址读取。捕获或者跟踪环路得到的码相位估计整数部分作为C/A码查找表的地址进行寻址，不足一个码片部分通过数据通道控制字TapOffset控制左/右移动数据实现。码相位计数器作为增计数器，由高速时钟驱动，当计数器计到1023时，生成码周期脉冲，用于跟踪电路其他模块的触发信号，如20ms计数器、整数ms相干累加等，同时重置码相位计数器为0。真实的产品中又是如何完成的呢？是通过移位方式还是内插方式产生本地适配的伪码呢？请在本系列文章中寻找答案。

2025年3月，再次续写本文，历经三年多的准备，这次将给出所有跟踪模块涉及的代码！后续文章将记录一边开发产品，一边展示代码并进行讲解的过程。本系列文章一出，信号跟踪算法明矣！本人特地为内容录制了多集视频，均放在视频号算法工匠中！

本谈将以案例展示的方式进行代码的讲解，方便大家从工程角度进行学习和消化！本篇文字接近八千字，有很多代码内容，非专业人士误入！视频中所用的代码在文章的下半部分展示了！本人已经将采样文件上传至百度网盘，永久共享，方便大家研究使用。  

跟踪内容综述！

• 载波跟踪：精确跟踪信号的载波频率和相位，消除多普勒频偏。需要用锁相环实现。

• 码跟踪：精确跟踪 C/A 码的相位，确保码同步。需要码跟踪环实现。

• 导航数据解调：从跟踪后的信号中解调出导航数据。

锁相环（PLL）：

• 原理：通过相位检测器、环路滤波器和压控振荡器（VCO）实现载波频率和相位的精确跟踪。

• 实现步骤：

1. 相位检测器：计算接收信号与本地载波信号的相位误差。

2. 环路滤波器：对相位误差进行滤波，生成控制信号。

3. 压控振荡器（VCO）：根据控制信号调整本地载波频率。

延迟锁定环（DLL）：

• 原理：通过早迟相关器、环路滤波器和码生成器实现 C/A 码相位的精确跟踪。

• 实现步骤：

   1. 早迟相关器：生成早码、准时码和迟码，计算相关值。

   2. 环路滤波器：对相关值进行滤波，生成控制信号。

   3. 码生成器：根据控制信号调整本地 C/A 码相位。

导航数据解调：

对剥离后的信号进行积分和判决，提取导航数据，然后进行定位解算。

跟踪模块的公开代码展示！

很多年前，我们就在苦苦追赶！

这里的展示有助于了解国内外代码编写的差异！

中文注释为本人所加，助力大家读懂代码！

黄色字体需要引起重视！

function [trackResults, channel]= tracking(fid, channel, settings)

% Performs code and carrier tracking for all channels.

%[trackResults, channel] = tracking(fid, channel, settings)

%   Inputs:

%       fid             - file identifier of the signal record.

%       channel     - PRN, carrier frequencies and code phases of all

%                       satellites to be tracked (prepared by preRum.m from

%                       acquisition results).

%       settings     - receiver settings.

%   Outputs:

%       trackResults    - tracking results (structure array). Contains

%                       in-phase prompt outputs and absolute spreading

%                       code's starting positions, together with other

%                       observation data from the tracking loops. All are

%                       saved every millisecond.

%-------------------------------------

%            SoftGNSS v3.0

% Copyright (C) Dennis M. Akos

% Written by Darius Plausinaitis and Dennis M. Akos

% Based on code by DMAkos Oct-1999

%--------------------------------------

%This program is free software。

%you can redistribute it and/or

%modify it under the terms of the GNU General Public License

%as published by the Free Software Foundation; either version 2

%of the License, or (at your option) any later version.

%

%This program is distributed in the hope that it will be useful,

%but WITHOUT ANY WARRANTY; without even the implied warranty of

%MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

%GNU General Public License for more details.

%You should have received a copy of the GNU General Public License

%along with this program; if not, write to the Free Software

%Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301,

%USA.

%CVS record:

%$Id: tracking.m,v 1.14.2.31 2006/08/14 11:38:22 dpl Exp $

%% Initialize result structure =====

% Channel status

trackResults.status   = '-';      

% No tracked signal, or lost lock

% The absolute sample in the record of the C/A code start:

trackResults.absoluteSample = zeros(1, settings.msToProcess);

% Freq of the C/A code:

trackResults.codeFreq       = inf(1, settings.msToProcess);

% Frequency of the tracked carrier wave:

trackResults.carrFreq       = inf(1, settings.msToProcess);

% Outputs from the correlators (In-phase):

trackResults.I_P            = zeros(1, settings.msToProcess);

trackResults.I_E            = zeros(1, settings.msToProcess);

trackResults.I_L            = zeros(1, settings.msToProcess);

% Outputs from the correlators (Quadrature-phase):

trackResults.Q_E            = zeros(1, settings.msToProcess);

trackResults.Q_P            = zeros(1, settings.msToProcess);

trackResults.Q_L            = zeros(1, settings.msToProcess);

% Loop discriminators

trackResults.dllDiscr       = inf(1, settings.msToProcess);

trackResults.dllDiscrFilt   = inf(1, settings.msToProcess);

trackResults.pllDiscr       = inf(1, settings.msToProcess);

trackResults.pllDiscrFilt   = inf(1, settings.msToProcess);

%--- Copy initial settings for all channels ------------

trackResults = repmat(trackResults, 1, settings.numberOfChannels);         %%将trackResults复制（分配）到各个Channels。

%% Initialize tracking variables ======

codePeriods = settings.msToProcess;     

% For GPS one C/A code is one ms

%--- DLL variables --------

% Define early-late offset (in chips)

earlyLateSpc = settings.dllCorrelatorSpacing;                              

%%步长设置。默认为0.5chip

% Summation interval

PDIcode = 0.001;

% Calculate filter coefficient values %%计算滤波器参数for PLL

[tau1code, tau2code] = calcLoopCoef(settings.dllNoiseBandwidth, settings.dllDampingRatio, 1.0);

%--- PLL variables ---------

% Summation interval

PDIcarr = 0.001;

% Calculate filter coefficient values

%%计算滤波器参数for DLL

[tau1carr, tau2carr] = calcLoopCoef(settings.pllNoiseBandwidth, ...

                                    settings.pllDampingRatio, 0.25);

hwb = waitbar(0,'Tracking...');

%% Start processing channels ====r 

channelNr = 1:settings.numberOfChannels

    % Only process if PRN is non zero (acquisition was successful)

    if (channel(channelNr).PRN ~= 0)

%%仅对捕获到了卫星进行跟踪

        % Save additional information - each channel's tracked PRN

        trackResults(channelNr).PRN     = channel(channelNr).PRN; %%分配卫星号

        % 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). In addition skip through that data file to start at the

        % appropriate sample (corresponding to code phase). Assumes sample

        % type is schar (or 1 byte per sample) 

fseek(fid, settings.skipNumberOfBytes + channel(channelNr).codePhase-1,  'bof');

%%寻找要处理的采样数据文件起始位置

        % Get a vector with the C/A code sampled 1x/chip

        caCode = generateCAcode(channel(channelNr).PRN); %%产生该通道对应的CA码

        % Then make it possible to do early and late versions

        caCode = [caCode(1023) caCode caCode(1)];

%% 拓展CA码，以方便产生超前滞后两路信号

% 这个操作很实用！

        %--- Perform various initializations ---

        % define initial code frequency basis of NCO

        codeFreq      = settings.codeFreqBasis; %%default -> 1.023MHz

        % define residual code phase (in chips)

        remCodePhase  = 0.0;

        % define carrier frequency which is used over whole tracking period

        carrFreq      = channel(channelNr).acquiredFreq; %%读取捕获阶段得到的频率值

        carrFreqBasis = channel(channelNr).acquiredFreq;

        % define residual carrier phase

        remCarrPhase  = 0.0;

        %code tracking loop parameters

        oldCodeNco   = 0.0; %%DLL环路参数

        oldCodeError = 0.0;

        %carrier/Costas loop parameters

        oldCarrNco   = 0.0; %%PLL环路参数

        oldCarrError = 0.0;

%=== Process the number of specified code periods  ===

        for loopCnt =  1:codePeriods %%处理的周期数,即处理几个ms的数据

%% GUI update ----------

%% GUI进度条显示

            % The GUI is updated every 50ms. This way Matlab GUI is still

            % responsive enough. At the same time Matlab is not occupied

            % all the time with GUI task.

            if (rem(loopCnt, 50) == 0)

                try

                    waitbar(loopCnt/codePeriods, ...

                            hwb, ...

                            ['Tracking: Ch ', int2str(channelNr), ...

                            ' of ', int2str(settings.numberOfChannels), ...

                            '; PRN#', int2str(channel(channelNr).PRN), ...

                            '; Completed ',int2str(loopCnt), ...

                            ' of ', int2str(codePeriods), ' msec']);                       

                catch

                    % The progress bar was closed. It is used as a signal

                    % to stop, "cancel" processing. Exit.

                    disp('Progress bar closed, exiting...');

                    return

                end

            end

%% Read next block of data -------

            % Find the size of a "block" or code period in whole samples

            % Update the phasestep based on code freq (variable) and

            % sampling frequency (fixed)

            codePhaseStep = codeFreq / settings.samplingFreq;              

blksize = ceil((settings.codeLength-remCodePhase) / codePhaseStep);

%%%%%%%%%%%%

            % Read in the appropriate number of samples to process this

            % interation 

            [rawSignal, samplesRead] = fread(fid,  blksize, settings.dataType);

            rawSignal = rawSignal';  %transpose vector

            % If did not read in enough samples, then could be out of 

            % data - better exit 

            if (samplesRead ~= blksize)

disp('Not able to read the specified number of samples  for tracking, exiting!')

                fclose(fid);

                return

            end

%% Set up all the code phase tracking information ---   %%%%%%%%%%%%%%%%%%%%%%%%%%

            % Define index into early code vector

            tcode       = (remCodePhase-earlyLateSpc) : ...

                          codePhaseStep : ...

 ((blksize-1)*codePhaseStep+remCodePhase-earlyLateSpc);

            %% 产生超前延迟及本地三路的CA码            

            tcode2      = ceil(tcode) + 1;

            earlyCode   = caCode(tcode2);

            % Define index into late code vector

            tcode       = (remCodePhase+earlyLateSpc) : ...

                          codePhaseStep : ...

 ((blksize-1)*codePhaseStep+remCodePhase+earlyLateSpc);

            tcode2      = ceil(tcode) + 1;

            lateCode    = caCode(tcode2);

            % Define index into prompt code vector

            tcode       = remCodePhase :  codePhaseStep : ...

                          ((blksize-1)*codePhaseStep+remCodePhase);

            tcode2      = ceil(tcode) + 1;

            promptCode  = caCode(tcode2);

            remCodePhase = (tcode(blksize) + codePhaseStep) - 1023.0;

%% Generate the carrier frequency to mix the signal to baseband ---

            time    = (0:blksize) ./ settings.samplingFreq;

% Get the argument to sin/cos functions

            trigarg = ((carrFreq * 2.0 * pi) .* time) + remCarrPhase;

            remCarrPhase = rem(trigarg(blksize+1), (2 * pi));

% Finally compute the signal to mix the collected data to bandband

            carrCos = cos(trigarg(1:blksize));

%% 产生本地正余弦波

            carrSin = sin(trigarg(1:blksize));

%% Generate the six standard accumulated values --------

            % First mix to baseband

            qBasebandSignal = carrCos .* rawSignal;

            %% 混频

            iBasebandSignal = carrSin .* rawSignal;

            % Now get early, late, and prompt values for each 相关运算

            I_E = sum(earlyCode  .* iBasebandSignal);

            Q_E = sum(earlyCode  .* qBasebandSignal);

            I_P = sum(promptCode .* iBasebandSignal);

            Q_P = sum(promptCode .* qBasebandSignal);

            I_L = sum(lateCode   .* iBasebandSignal);

            Q_L = sum(lateCode   .* qBasebandSignal);

%% Find PLL error and update carrier NCO ----     

%% 鉴频鉴相

% Implement carrier loop discriminator (phase detector)

            carrError = atan(Q_P / I_P) / (2.0 * pi);

% 这里的鉴相公式有个除法运算哦，很实用！书本上可不会这么写！

% Implement carrier loop filter and generate NCO command

            carrNco = oldCarrNco + (tau2carr/tau1carr) * ...

                (carrError - oldCarrError) + carrError * (PDIcarr/tau1carr);

            oldCarrNco   = carrNco;

            oldCarrError = carrError;

            % Modify carrier freq based on NCO command

            carrFreq = carrFreqBasis + carrNco;

            trackResults(channelNr).carrFreq(loopCnt) = carrFreq;

%% Find DLL error and update code NCO -------

            codeError = (sqrt(I_E * I_E + Q_E * Q_E) - sqrt(I_L * I_L + Q_L * Q_L)) / ...

                (sqrt(I_E * I_E + Q_E * Q_E) + sqrt(I_L * I_L + Q_L * Q_L));

            % Implement code loop filter and generate NCO command

            codeNco = oldCodeNco + (tau2code/tau1code) * ...

                (codeError - oldCodeError) + codeError * (PDIcode/tau1code);

            oldCodeNco   = codeNco;

            oldCodeError = codeError;

            % Modify code freq based on NCO command

            codeFreq = settings.codeFreqBasis - codeNco;

            trackResults(channelNr).codeFreq(loopCnt) = codeFreq;

%% Record various measures to show in postprocessing -----

            % Record sample number (based on 8bit samples)

            trackResults(channelNr).absoluteSample(loopCnt) = ftell(fid);

            trackResults(channelNr).dllDiscr(loopCnt)       = codeError;

            trackResults(channelNr).dllDiscrFilt(loopCnt)   = codeNco;

            trackResults(channelNr).pllDiscr(loopCnt)       = carrError;

            trackResults(channelNr).pllDiscrFilt(loopCnt)   = carrNco;

            trackResults(channelNr).I_E(loopCnt) = I_E;

            trackResults(channelNr).I_P(loopCnt) = I_P;

            trackResults(channelNr).I_L(loopCnt) = I_L;

            trackResults(channelNr).Q_E(loopCnt) = Q_E;

            trackResults(channelNr).Q_P(loopCnt) = Q_P;

            trackResults(channelNr).Q_L(loopCnt) = Q_L;

        end % for loopCnt

        % If we got so far, this means that the tracking was successful

        % Now we only copy status, but it can be update by a lock detector

        % if implemented

        trackResults(channelNr).status  = channel(channelNr).status;     

    end % if a PRN is assigned

end % for channelNr 

% Close the waitbar

close(hwb)

看看本人二十多年前写的跟踪程序吧！

需要了解GPS信号采样文件的参数和格式！

在捕获程序中有对应的参数设置！

主要是采样率和中频！