GPS C/A码捕获与跟踪的MATLAB实现：从Gold码原理到环路参数调试

简介：一套完整的GPS信号处理MATLAB工程，覆盖C/A码生成、GPS信号模拟、卫星信号捕获与跟踪全套流程，运行后可输出捕获结果、相位起始位置、跟踪偏差及其震荡曲线。工程采用模块化源码组织，基于FFT相关峰值搜索与早迟门跟踪环路，各功能独立封装便于单独调试与参数修改，适合通信导航专业学生、研究者以及GPS接收机开发工程师用于教学演示或二次开发。压缩包共11个文件，其中7个.m脚本分别实现主程序、捕获、跟踪、码表生成、环路系数计算等功能，另有4个.asv自动备份文件，整体仅15KB，轻量易用。已有578人学习下载。通过调节噪声强度等参数，可模拟不同信噪比环境，直观观察相关峰定位、偏差收敛过程及震荡形态，从而深入理解C/A码特性与GPS信号处理各环节的算法实现，对学习GPS接收机原理和信号处理流程具有较高参考价值。

1. 用 MATLAB 跑通 GPS C/A 码捕获与跟踪：先解决为什么跟踪比捕获难

GPS 的 C/A 码捕获，很多教材会把它讲成一个二维搜索问题：在 1023 个码相位和几十个多普勒频点里找峰值。真正用 MATLAB 把它跑通之后你会发现，捕获只是拿到“相位起始位置和粗略载波频率”的入场券，后面的跟踪才是决定接收机能不能稳定输出的部分。捕获结果里有一个相位起始位置，跟踪环路里有一套偏差序列，这两段能不能衔接上，直接决定你是在做演示还是做接收机原型。

这份笔记按“产生码 → 捕获 → 跟踪 → 排错 → 验证”的顺序来写。适合两类人：一类是刚接触 GPS/GNSS 仿真的学生，需要一套能跑通的 MATLAB 主链路；另一类是已经把捕获跑通、但跟踪环一直发散的工程师，可以直接跳到第 5 章看边界条件和参数设置。所有代码都按常见做法写成最小可运行形式，你只需要把数据读取函数换成自己的文件路径。

2. 产生 C/A 码：从 Gold 码原理到一个可复现的 MATLAB 函数

2.1 为什么 C/A 码要用 Gold 码：从两个 G2 抽头说起

GPS L1 频段的 C/A 码是 Gold 码，码率 1.023 MHz，周期 1023 码片，正好 1 ms。它不是单个 m 序列，而是由两个 10 级移位寄存器 G1 和 G2 组合出来的。G2 的输出不是直接取最后一级，而是从两个抽头位置异或得到，不同卫星对应不同的抽头组合。这套设计让不同卫星码之间的互相关峰值有上界，不会因为互相关高而互相误触发。

实现时最容易被带偏的一点是：很多资料只讲“G1 和 G2 反馈多项式”，却不说 G2 抽头值的决定作用。同一个 G2 寄存器，选不同抽头异或，出来的序列相位完全不同，生成的 C/A 码也就不同。所以在写生成函数之前，先把 PRN 号对应的抽头表准备好，否则你生成的自相关看起来正常，和真实 GPS 数据却始终对不上。

具体细节也不复杂：G1 的反馈抽头通常用第 3 级和第 10 级；G2 的反馈抽头在第 2、3、6、8、9、10 级。初始状态全部置 1。每个码片输出是 G1 第 10 级和 G2 两个抽头异或值的异或。下面代码里给出的抽头对是 1 到 8 号卫星的常用值，完整 32 个抽头表可以在 GPS ICD 文档或开源接收机代码里找到。

2.2 一个可直接跑的 C/A 码生成函数（MATLAB）

下面这个函数完整可运行，返回 1023×1 的 ±1 序列。注意我刻意返回 ±1 而不是 0/1，是因为后续捕获的匹配滤波和相关峰检测需要双极性信号。实际工程中，把 0/1 映射成 ±1 并不会影响相关峰位置，但会直接影响相关输出的绝对幅度，所以统一在生成端做好。

function caCode = generateCaCode(prn) % generateCaCode 生成 GPS C/A 码 % 输入 prn：卫星号，示例支持 1~8 % 输出 caCode：1023x1 的 +/-1 序列 tapTable = [2 6; 3 7; 4 8; 5 9; 1 9; 2 10; 1 8; 2 9]; % PRN 1-8 if prn < 1 || prn > 8 error('示例代码只支持 PRN 1~8，其余卫星请补抽头表'); end tapA = tapTable(prn, 1); tapB = tapTable(prn, 2); g1 = ones(1, 10); g2 = ones(1, 10); caCode = zeros(1, 1023); for i = 1:1023 g1out = g1(10); g2out = xor(g2(tapA), g2(tapB)); caCode(i) = 1 - 2 * xor(g1out, g2out); % 映射成 +1/-1 fb1 = xor(g1(3), g1(10)); % G1 反馈 g1 = [fb1, g1(1:9)]; fb2 = xor(xor(xor(xor(xor(g2(2), g2(3)), g2(6)), g2(8)), g2(9)), g2(10)); g2 = [fb2, g2(1:9)]; end caCode = caCode(:); end

循环里先取当前输出，再更新寄存器，顺序不能反。G2 抽头组合tapA、tapB异或得到的值，等效于 G2 序列某个特定相位；再和 G1 末级异或就得到当前码片。如果你把寄存器存成logical类型，xor也同样工作，但为了后续直接做乘法相关，这里保持双精度数组更省事。

2.3 验证生成结果：自相关尖峰与互相关低电平

生成之后先别急着进捕获，用相关函数验证码是否正确。C/A 码自相关在零移位是 1023，其他移位接近 -1；互相关通常在 ±65 附近。如果不满足这个形态，多半是抽头表或反馈多项式写错了。

ca1 = generateCaCode(1); ca2 = generateCaCode(2); [autoCorr, lags] = xcorr(ca1, ca1, 1022); [crossCorr, ~] = xcorr(ca1, ca2, 1022); figure('Name', 'C/A码验证'); subplot(2, 1, 1); plot(lags, autoCorr); title('PRN1 自相关'); ylabel('幅度'); subplot(2, 1, 2); plot(lags, crossCorr); title('PRN1-PRN2 互相关'); ylabel('幅度');

这里为了图形直观用了线性xcorr。严格讲，GPS 信号中的 C/A 码是周期序列，捕获阶段用的是循环相关；但线性相关足以暴露抽头错误。你可以再用max(abs(crossCorr))检查互相关峰值是否远小于 1023，这个值越大，越说明抽头表选得有问题。

这一步也是后续捕获的“本地码基准”。基准错了，捕获和跟踪全部白做。我一般会把 PRN 1 和 PRN 2 的码序列打印成图，确认互相关确实接近噪声电平，再继续往下走。

3. 捕获信号：用 FFT 并行搜索码相位与多普勒频移

3.1 捕获的本质：在 1023 个码相位和若干个多普勒频点上找能量峰值

接收机拿到的是经过下变频、采样的数字信号。信号里除了噪声，还有一颗卫星的 C/A 码序列，但这个码的起始位置完全未知。卫星相对接收机的运动还会带来多普勒频率偏移，L1 频段上最大可达 ±5 kHz 左右。所以捕获就是在一个二维平面上找峰值：一维是多普勒频率，另一维是码相位。码相位候选不是 1023 个码片，而是按采样率展开后的采样点总数。例如采样率 5.115 MHz，1 ms 内有 5115 个采样点，相位候选就有 5115 个。

串行捕获的思路是每个候选码相位都做一次相关，复杂度接近 O(N²)，在 MATLAB 里跑起来非常慢。工程上几乎都用 FFT 做并行码相位搜索：对固定多普勒频率，一次 FFT 相关就算出了所有码相位的相关值。多普勒维度仍然要扫描，但可以减小步进或者用粗频点减少数量。积分时间用 1 ms，对应一个完整 C/A 码周期，频率分辨率约 1 kHz，所以多普勒搜索步进在实际工程里取 500 Hz 或 250 Hz 都能用。

3.2 用 FFT 并行搜索码相位：MATLAB 实现

下面函数实现了一次 1 ms 数据上的并行码相位搜索。输入假定为基带复数信号，也就是已经去掉了中频载波；如果你输入的是真实前端输出的中频信号，需要先乘一个本地中频载波再进这个函数。函数返回最大相关幅度、对应的码相位起始采样点、以及多普勒频率。

function [peakVal, codePhaseSample, dopplerFreq] = acquireCa(signal, fs, prn, fdRange, fdStep) % signal : 1ms 基带复数采样，长度约 fs/1000 % fs : 采样率，单位 Hz % prn : 需要捕获的卫星号 % fdRange: 多普勒搜索范围，例如 [-5000 5000] % fdStep : 多普勒搜索步进，例如 500 或 250 % 输出 codePhaseSample 为 1ms 内的起始采样点（从 1 开始） N = length(signal); t = (0:N-1) / fs; caCode = generateCaCode(prn); fdList = fdRange(1):fdStep:fdRange(2); peakVal = 0; codePhaseSample = 1; dopplerFreq = 0; for k = 1:length(fdList) fd = fdList(k); carrier = exp(-1j * 2 * pi * fd * t); y = signal .* carrier; % 将 1023 码片展开到采样域 codeIndex = mod(floor(t * 1.023e6), 1023) + 1; localCode = caCode(codeIndex); % 循环互相关 corr = ifft(fft(y) .* conj(fft(localCode))); [p, idx] = max(abs(corr)); if p > peakVal peakVal = p; codePhaseSample = idx; dopplerFreq = fd; end end end

代码最核心的是corr = ifft(fft(y) .* conj(fft(localCode)))。fft(y)把输入信号变到频域，conj(fft(localCode))相当于本地码的匹配滤波，再做逆变换得到所有循环移位的相关值。循环相关的含义是：码相位可以在 1 ms 内任意循环，正好匹配 C/A 码周期。如果本地码不是周期延拓而是线性序列，峰值位置会丢失大半。

参数上，fdStep选择 500 Hz 时，1 ms 相干积分的频率响应主瓣宽度约 1 kHz，500 Hz 步进不会漏掉主瓣峰。如果想捕获弱信号或对频率精度要求高，改成 250 Hz，但计算量翻倍。fs不建议低于 2.046 MHz，否则一个码片不足两个采样点，相关器分辨率不够。实际接收机前端常用 5.115 MHz，对应每码片 5 个采样点，这个函数直接用这个采样率即可。

3.3 输出捕获结果：峰位图、相位起始位置与多普勒频率

调用上面的函数后，你需要做两件事：门限判决和出图。门限不能只看峰值绝对值，要与噪声基底比。一般做法是把当前相关功率与相关输出的平均功率作比，超过 8~12 dB 判为捕获成功。若使用多帧非相干累加，门限可以进一步降低。

% 读取 1ms 数据，readIfData 是自实现的读取函数，返回复数采样 signal = readIfData('my_gps_data.bin', fs); % 捕获卫星 1 [peakVal, codePhaseSample, dopplerFreq] = acquireCa(signal, fs, 1, [-5000 5000], 500); % 门限判断：与相关谱平均水平比较 N = length(signal); t = (0:N-1)/fs; caCode = generateCaCode(1); localCode = caCode(mod(floor(t*1.023e6), 1023)+1); corr = ifft(fft(signal .* exp(-1j*2*pi*dopplerFreq*t)) .* conj(fft(localCode))); noisePower = mean(abs(corr).^2); peakPower = max(abs(corr).^2); fprintf('峰均比 = %.1f dB\n', 10*log10(peakPower/noisePower)); % 换算成码片相位 codePhaseChip = (codePhaseSample - 1) * 1.023e6 / fs; fprintf('捕获结果：码相位起始 = %.2f 码片，多普勒 = %.2f Hz\n', codePhaseChip, dopplerFreq); plot(abs(corr)); hold on; plot(codePhaseSample, sqrt(peakPower), 'ro'); xlabel('1ms 内采样点'); ylabel('相关幅度');

codePhaseSample是采样点索引，如果你后续要送给跟踪环，直接用这个采样点即可；如果想看码相位，用(idx-1)*1.023e6/fs换算成码片。注意这里输出的是“一个码周期内的起始位置”，不包含整周模糊度；C/A 码每 1 ms 重复一次，所以这个值是 0~1023 码片范围内的小数。

门限的选取和信号强度强相关。我的经验是：仿真纯净信号下峰均比可以到 20 dB 以上，而真实天线数据经常只有 8~10 dB，所以不要拿仿真门限去接真实数据，否则虚警会高得离谱。这一步的相位起始位置，就是标题里要求输出的第一个关键量。

4. 跟踪信号：码环与载波环联动，输出误差序列与偏差震荡

4.1 跟踪环路的基本结构：码环和载波环如何协同

捕获给出的是粗估计，误差通常在几十 Hz 和一个码片以内。跟踪阶段要把这两个残差压到接近零，并且持续输出。跟踪由两个环组成：载波环负责把残余多普勒和载波相位剥离干净，码环负责让本地码相位对齐输入码相位。码环通常用超前、即时、滞后三路相关器，取超前减滞后功率差作为误差；载波环用 Costas 环，取atan(QP/IP)作为相位误差，可以避免 180° 相位反转带来的比特翻转问题。

需要注意，码环和载波环不是独立的。载波环输出的频率估计会通过多普勒牵引影响码速率：码率约等于 1.023 MHz 乘上频率偏差与 L1 频率之比。如果忽略这一项，动态场景下码环会持续偏出去，导致跟踪偏差震荡整体抬升。这也是新手从“捕获成功”到“跟踪发散”之间最常见的坑。

4.2 一个 1ms 更新一次的跟踪循环（MATLAB）

下面的trackCa函数把整个环路更新写在一个循环里，每个 1 ms 处理一块信号，输出每一毫秒的码环鉴别器和载波环鉴别器。环路滤波器用二阶积分结构，参数先用带宽和阻尼系数归一化。实际使用时需要根据输入信号的动态和载噪比微调。

function [carrierErr, codeErr] = trackCa(signal, fs, prn, initCodePhase, initDoppler) % 输入信号为连续多毫秒基带复数采样 % initCodePhase: 捕获得到的起始采样点 % initDoppler : 捕获得到的多普勒频率初值 msLen = round(fs/1000); Nms = floor(length(signal)/msLen); ca = generateCaCode(prn); zeta = 0.7; bwC = 1; % 码环带宽，单位 Hz bwP = 15; % 载波环带宽，单位 Hz T = 0.001; % 环路更新周期 wnC = bwC*4/(1+4*zeta^2); wnP = bwP*4/(1+4*zeta^2); alphaC = 2*zeta*wnC*T; betaC = (wnC*T)^2; alphaP = 2*zeta*wnP*T; betaP = (wnP*T)^2; codePhase = initCodePhase; fd = initDoppler; carrierPhase = 0; codeFreqAcc = 0; freqAcc = 0; dSamp = 0.5 * fs / 1.023e6; carrierErr = zeros(Nms-1, 1); codeErr = zeros(Nms-1, 1); for k = 1:Nms seg = signal((k-1)*msLen+1 : k*msLen); t = (0:msLen-1)/fs; carrier = exp(-1j*(2*pi*fd*t + carrierPhase)); x = seg .* carrier; % 码率随载波多普勒牵引 codeNcoStep = fs / (1.023e6 * (1 + fd/1.57542e9)); ixP = mod(floor(codePhase + (0:msLen-1)*codeNcoStep), 1023) + 1; ixE = mod(floor(codePhase - dSamp + (0:msLen-1)*codeNcoStep), 1023) + 1; ixL = mod(floor(codePhase + dSamp + (0:msLen-1)*codeNcoStep), 1023) + 1; lcP = ca(ixP); lcE = ca(ixE); lcL = ca(ixL); IP = sum(real(x).*lcP); QP = sum(imag(x).*lcP); IE = sum(real(x).*lcE); QE = sum(imag(x).*lcE); IL = sum(real(x).*lcL); QL = sum(imag(x).*lcL); codeDisc = (IE^2+QE^2 - IL^2-QL^2) / (IE^2+QE^2 + IL^2+QL^2 + eps); carrierDisc = atan(QP/IP); codeFreqAcc = codeFreqAcc + betaC*codeDisc; codePhase = codePhase + codeNcoStep*msLen + alphaC*codeDisc + codeFreqAcc; freqAcc = freqAcc + betaP*carrierDisc; fd = fd + alphaP*carrierDisc + freqAcc; carrierPhase = mod(carrierPhase + 2*pi*fd*msLen/fs, 2*pi); if k > 1 carrierErr(k-1) = carrierDisc; codeErr(k-1) = codeDisc; end end end

这段代码是工程里常见的离散化结构。注意环路滤波器的系数都乘以了T = 0.001，这是最容易写错的一步；少了T，环路增益会大几百倍，跟踪必然发散。codeNcoStep用fd/1.57542e9修正码率，是因为 GPS L1 载波频率是 1575.42 MHz，码率和载波频率受同一多普勒影响，忽略这个修正项，动态场景下码环会一直有残差。

参数说明：bwC=1 Hz、bwP=15 Hz是静态场景的常用起点。带宽越大，动态响应越快，但噪声越差；带宽越小，输出越干净，但跟不上加加速度。dSamp=0.5是半码片间距，弱信号下可以放宽到 1 码片。这些参数就是决定“偏差震荡”大小的关键，后面的第 5 章会专门讲它们带来的坑。

4.3 输出跟踪偏差和偏差震荡：记录序列并统计

跟踪完成之后，carrierErr和codeErr就是你要的跟踪偏差。偏差震荡通常用峰峰值或标准差描述。下面代码直接出图并打印两个指标。

[carrierErr, codeErr] = trackCa(signal, fs, 1, codePhaseSample, dopplerFreq); figure('Name', '跟踪误差输出'); subplot(2, 1, 1); plot(codeErr, '.-'); ylabel('码相位误差(码片)'); grid on; subplot(2, 1, 2); plot(carrierErr, '.-'); ylabel('载波相位误差(rad)'); grid on; fprintf('码偏差震荡 std = %.4f 码片，峰峰值 = %.4f 码片\n', ... std(codeErr), max(codeErr)-min(codeErr)); fprintf('载波偏差震荡 std = %.4f rad，峰峰值 = %.4f rad\n', ... std(carrierErr), max(carrierErr)-min(carrierErr));

若输入是高斯白噪声环境下的合成信号，跟踪稳定后codeErr会围绕 0 振荡，std一般小于 0.05 码片；载波相位误差的std应该小于 0.1 rad 量级。如果你看到误差序列不断漂移或周期性跳变，先检查是否把捕获的初值传给了跟踪，再检查环路带宽是不是太小。

这里统计的就是标题里说的“偏差震荡”。好的跟踪输出，应该是偏差均值接近 0、震荡幅度随带宽变窄而变小。注意不要直接把鉴别器输出当最终测量值，它只是闭环误差信号；真正用于定位的伪距，还要从码 NCO 相位和载波相位累加器里读。

5. 避坑：捕获跟踪链路中 5 个最常翻车的点

5.1 现象：捕获相关峰是平的，找不到尖峰

跑acquireCa之后，max(abs(corr))始终和均值差不多，二维搜索里没有明显凸起。最常见原因是信号端和本地码的码片速率没对齐，比如采样率 5.115 MHz，但本地码生成时用了 1.023 MHz 直接索引；另一种可能是输入信号实际是中频信号，没有做下变频，FFT 相关在搜索多普勒时没法抵消中频分量，相关结果自然被压平。

解决方式也很直接：先用一个已知的合成信号自测。生成 1 ms 纯 C/A 码信号，叠加一个 1 kHz 多普勒，再让捕获函数去恢复。如果合成信号能捕到，真实数据捕不到，就往数据格式、前端增益和中频下变频上查。不要一开始就怀疑算法，这一步能把问题快速框定在数据链路还是算法链路。

5.2 现象：捕获到的相位起始位置总是差半个码片

捕获峰值很高，但送进跟踪后收敛到邻近码片，或者和已知延迟对比时总是差 0.5 码片左右。原因是 FFT 循环相关的分辨率取决于采样率，采样点会落在码片边界上，峰值位置被量化到离真实相位最近的整采样点；另一个原因是本地码索引用了floor，相位本身就只能表示到采样点。

解决方法是提高采样率，或者在峰值附近做抛物线插值。codePhaseSample是采样点索引，无法表达小数码片；跟踪初值可以接受这个误差，但如果你要做高精度测距，必须补插值。在峰值点和左右相邻点之间做抛物线插值，可以让相位精度从 1 个采样点提升到 0.1 个采样点以内。

5.3 现象：跟踪环路在强信号下还是发散

信号很强，相关峰也很尖，但trackCa跑到几十毫秒就开始震荡甚至失锁。我见过的最多原因是环路滤波器系数单位弄错了。很多资料把带宽直接代进系数，却忘了离散化后系数与更新周期T强相关。alpha和beta必须与T配套；如果写成alpha = 2*zeta*wn少了乘T，环路增益大了几百倍，发散是必然的。

排查时先检查alpha、beta的计算是否包含T = 0.001。如果还是发散，把载波环带宽从 15 Hz 降到 5 Hz 试试。强信号下发散基本不是噪声问题，而是相位裕度或增益问题。

5.4 现象：偏差震荡很大，均值却正常

codeErr在 0 附近，但std很大，超过 0.2 码片；carrierErr也在 0 附近但抖动达到几十度。这可能不是环路本身的问题，而是输入信号的载噪比太低，或者相干积分长度不是 1 ms。在真实弱信号场景下，误差抖动和环路带宽直接相关：带宽越大，输出越脏。

解决思路是让带宽跟着信号强度走。常见做法是先用 15 Hz 载波环拉入跟踪，稳定后切到 2 Hz 窄带宽保持精度。这个“先宽后窄”的技巧能显著降低偏差震荡，也是接收机里自动增益控制和载噪比估计之后最常用的策略。

5.5 现象：FFT 相关出现两个几乎一样的峰

输出相关谱里有相邻的两个峰，高度接近，导致max返回的索引来回跳，造成相位起始位置抖动。当码相位差接近整数码片，而且采样率刚好等于码速率的整数倍时，码片边界的采样点会跨越两个码片，产生栅瓣。另一个原因是多普勒频率没有对准，残余频率使相关峰展宽，看起来像两个峰。

解决方法是把fdStep改小到 250 Hz，或者在峰值判决前对相关幅度做三采样点抛物线插值，再取最大位置。频率维度也可以做插值，用相邻两个频点的幅度比来估计真实多普勒，这样捕获输出的频率更准，跟踪阶段收敛更快。

6. 用仿真与真实数据双向验证：把捕获参数和跟踪参数固化下来

跟踪环调通后，不要急着声称算法完成。我的验证顺序是先“自产自销”：在 MATLAB 里生成一个已知码相位和多普勒的 C/A 码信号，加入高斯噪声，用捕获和跟踪去恢复，然后比较误差。这一步能暴露九成以上的逻辑错误。

真实数据可以从常见的公开 GPS 中频数据集里找，采样率通常落在 5 MHz 到 20 MHz 之间。注意先确认数据文件是复数还是实数、字节序、是否有直流偏置。很多“捕不到”的问题其实是读错数据格式，和算法没关系。

一个非常值得固化的细节是：把捕获阶段的codePhaseSample和dopplerFreq直接作为跟踪初值。这个道理看起来显而易见，但很容易被忽略。我见过太多人把初始化写死成 0，导致跟踪环花了很久去追，甚至追不上。另外，环路带宽最好动态切换：先用 15 Hz 拉入，500 ms 之后切成 2 Hz，可以明显减小第 4.3 节统计的偏差震荡。切换点可以用abs(carrierDisc)连续小于 0.1 rad 作为锁定条件，这样比固定时间点更可靠。

我自己有一次把一个项目里的 C/A 码捕获做好之后，给跟踪写了一个零初值初始化。当时用合成信号发现跟踪能收敛，只是慢，就没细想；换到真实数据后，跟踪一直在 180° 附近转，查了大半天才意识到是初值相位模糊。后来我养成了习惯：所有调试脚本都从捕获输出里取初值，并且把初值来源打印出来，省下无数“玄学”排错时间。希望这个习惯也能帮到你。

本文还有配套的精品资源，点击获取